Privratsky-20210215 - Sub cluster the “Macro” population
Hélène Fradin Kirshner
04/13/2021
source conda/bin/activate privratsky rmarkdown::render(‘./9_Subcluster_Macrophages/9_Subcluster_Macrophages.Rmd’)
Changes in myeloid and kidney cells after CLP - Analysis of 2 x 10X scRNA-seq samples from 2 pools of WT mice (3 Sham + 3 CLP): comparison of gene expression in different cell populations
Upregulation of anti-inflammatory mediators in macrophage population (SOCS3, IL-1R2, IL1rn) and downregulation of pro-inflammatory genes
indir <- "./processedData/7_DEA_C24_vs_C0"
outdir <- "./processedData/9_Subcluster_Macrophages"
dir.create(outdir, recursive = T)library(Seurat)
library(ggplot2)
annotated <- readRDS(paste0(indir, "/3.annotated.rds"))
annotated## An object of class Seurat
## 24399 features across 18055 samples within 2 assays
## Active assay: RNA (22399 features, 0 variable features)
## 1 other assay present: integrated
## 2 dimensional reductions calculated: pca, umapSubset to Macrophages and sub-cluster
annotated## An object of class Seurat
## 24399 features across 18055 samples within 2 assays
## Active assay: RNA (22399 features, 0 variable features)
## 1 other assay present: integrated
## 2 dimensional reductions calculated: pca, umapIdents(annotated) <- "annotation.1"
macrophages <- subset(annotated, idents = "Macro")
macrophages## An object of class Seurat
## 24399 features across 588 samples within 2 assays
## Active assay: RNA (22399 features, 0 variable features)
## 1 other assay present: integrated
## 2 dimensional reductions calculated: pca, umaptable(macrophages$annotation.1)##
## Endo Podo PT-S1 PT-S2 PT-S3
## 0 0 0 0 0
## LOH Distal Conv T Conn Tubule CD-PC CD-IC
## 0 0 0 0 0
## CD-Trans CD-IM Fib Macro PMN
## 0 0 0 588 0
## B lymph NK
## 0 0DefaultAssay(macrophages) <- "RNA"
macrophages.list <- SplitObject(macrophages, split.by = "sample.id")
macrophages.list## $C0
## An object of class Seurat
## 24399 features across 267 samples within 2 assays
## Active assay: RNA (22399 features, 0 variable features)
## 1 other assay present: integrated
## 2 dimensional reductions calculated: pca, umap
##
## $C24
## An object of class Seurat
## 24399 features across 321 samples within 2 assays
## Active assay: RNA (22399 features, 0 variable features)
## 1 other assay present: integrated
## 2 dimensional reductions calculated: pca, umapfor (i in 1:length(macrophages.list)) {
macrophages.list[[i]] <- FindVariableFeatures(macrophages.list[[i]],
selection.method = "vst", nfeatures = 2000, verbose = FALSE)
}k.filter <- min(200, min(sapply(macrophages.list, ncol)))
k.filter## [1] 200macrophages.anchors <- FindIntegrationAnchors(object.list = macrophages.list,
dims = 1:30, k.filter = k.filter)macrophages.integrated <- IntegrateData(anchorset = macrophages.anchors,
dims = 1:30)DefaultAssay(macrophages.integrated) <- "integrated"
macrophages.integrated <- ScaleData(macrophages.integrated, verbose = T)dim(macrophages.integrated)## [1] 2000 588macrophages.integrated <- RunPCA(macrophages.integrated, npcs = 50,
verbose = T)system.time(macrophages.integrated <- JackStraw(macrophages.integrated,
num.replicate = 100, dims = 50))## user system elapsed
## 74.031 0.002 74.104macrophages.integrated <- ScoreJackStraw(macrophages.integrated,
dims = 1:50)
j <- JackStrawPlot(macrophages.integrated, dims = 1:50)
pdf(paste0(outdir, "/1_JackStrawPlot_macrophages.pdf"))
j
dev.off()## png
## 2macrophages.integrated <- FindNeighbors(macrophages.integrated,
dims = 1:50)
macrophages.integrated <- RunUMAP(macrophages.integrated, dims = 1:50)for (i in seq(0, 2, 0.1)) {
macrophages.integrated <- FindClusters(macrophages.integrated,
resolution = i, verbose = F)
macrophages.integrated[[paste0("macrophages_subclusters_res.",
i)]] <- macrophages.integrated$seurat_clusters
}head(macrophages.integrated[[]])## orig.ident nCount_RNA nFeature_RNA percent.mito sample.id
## AAAGGGCAGTCACGAG--C0 C0 1233 661 8.191403 C0
## AAAGTCCCAGGTTCAT--C0 C0 3888 1562 9.027778 C0
## AAAGTCCTCCGAGGCT--C0 C0 4237 1508 2.336559 C0
## AACAAAGCACCTGCGA--C0 C0 3528 1579 4.365079 C0
## AACCTGACATGCGGTC--C0 C0 884 420 19.230769 C0
## AACGTCAAGGTTCATC--C0 C0 4528 1895 5.631625 C0
## integrated_snn_res.0 seurat_clusters
## AAAGGGCAGTCACGAG--C0 0 6
## AAAGTCCCAGGTTCAT--C0 0 2
## AAAGTCCTCCGAGGCT--C0 0 2
## AACAAAGCACCTGCGA--C0 0 2
## AACCTGACATGCGGTC--C0 0 3
## AACGTCAAGGTTCATC--C0 0 7
## integrated_snn_res.0.1 integrated_snn_res.0.2
## AAAGGGCAGTCACGAG--C0 0 0
## AAAGTCCCAGGTTCAT--C0 0 0
## AAAGTCCTCCGAGGCT--C0 0 0
## AACAAAGCACCTGCGA--C0 0 0
## AACCTGACATGCGGTC--C0 0 0
## AACGTCAAGGTTCATC--C0 1 2
## integrated_snn_res.0.3 integrated_snn_res.0.4
## AAAGGGCAGTCACGAG--C0 0 0
## AAAGTCCCAGGTTCAT--C0 0 0
## AAAGTCCTCCGAGGCT--C0 0 0
## AACAAAGCACCTGCGA--C0 0 0
## AACCTGACATGCGGTC--C0 0 0
## AACGTCAAGGTTCATC--C0 2 2
## integrated_snn_res.0.5 integrated_snn_res.0.6
## AAAGGGCAGTCACGAG--C0 0 0
## AAAGTCCCAGGTTCAT--C0 0 0
## AAAGTCCTCCGAGGCT--C0 0 0
## AACAAAGCACCTGCGA--C0 0 0
## AACCTGACATGCGGTC--C0 0 0
## AACGTCAAGGTTCATC--C0 2 2
## integrated_snn_res.0.7 integrated_snn_res.0.8
## AAAGGGCAGTCACGAG--C0 0 0
## AAAGTCCCAGGTTCAT--C0 0 0
## AAAGTCCTCCGAGGCT--C0 0 0
## AACAAAGCACCTGCGA--C0 0 0
## AACCTGACATGCGGTC--C0 0 0
## AACGTCAAGGTTCATC--C0 2 2
## integrated_snn_res.0.9 integrated_snn_res.1
## AAAGGGCAGTCACGAG--C0 0 0
## AAAGTCCCAGGTTCAT--C0 0 0
## AAAGTCCTCCGAGGCT--C0 0 0
## AACAAAGCACCTGCGA--C0 0 0
## AACCTGACATGCGGTC--C0 0 0
## AACGTCAAGGTTCATC--C0 2 2
## integrated_snn_res.1.1 integrated_snn_res.1.2
## AAAGGGCAGTCACGAG--C0 0 3
## AAAGTCCCAGGTTCAT--C0 0 0
## AAAGTCCTCCGAGGCT--C0 0 0
## AACAAAGCACCTGCGA--C0 0 0
## AACCTGACATGCGGTC--C0 0 0
## AACGTCAAGGTTCATC--C0 2 2
## integrated_snn_res.1.3 integrated_snn_res.1.4
## AAAGGGCAGTCACGAG--C0 2 2
## AAAGTCCCAGGTTCAT--C0 1 1
## AAAGTCCTCCGAGGCT--C0 1 1
## AACAAAGCACCTGCGA--C0 2 2
## AACCTGACATGCGGTC--C0 2 2
## AACGTCAAGGTTCATC--C0 3 3
## integrated_snn_res.1.5 integrated_snn_res.1.6
## AAAGGGCAGTCACGAG--C0 2 3
## AAAGTCCCAGGTTCAT--C0 0 0
## AAAGTCCTCCGAGGCT--C0 0 4
## AACAAAGCACCTGCGA--C0 2 4
## AACCTGACATGCGGTC--C0 2 3
## AACGTCAAGGTTCATC--C0 3 2
## integrated_snn_res.1.7 integrated_snn_res.1.8
## AAAGGGCAGTCACGAG--C0 3 3
## AAAGTCCCAGGTTCAT--C0 0 2
## AAAGTCCTCCGAGGCT--C0 4 2
## AACAAAGCACCTGCGA--C0 4 2
## AACCTGACATGCGGTC--C0 3 3
## AACGTCAAGGTTCATC--C0 1 7
## integrated_snn_res.1.9 integrated_snn_res.2
## AAAGGGCAGTCACGAG--C0 3 6
## AAAGTCCCAGGTTCAT--C0 2 2
## AAAGTCCTCCGAGGCT--C0 2 2
## AACAAAGCACCTGCGA--C0 2 2
## AACCTGACATGCGGTC--C0 8 3
## AACGTCAAGGTTCATC--C0 7 7
## integrated_snn_res.2.1 integrated_snn_res.2.2
## AAAGGGCAGTCACGAG--C0 26 29
## AAAGTCCCAGGTTCAT--C0 18 22
## AAAGTCCTCCGAGGCT--C0 18 22
## AACAAAGCACCTGCGA--C0 18 22
## AACCTGACATGCGGTC--C0 26 29
## AACGTCAAGGTTCATC--C0 36 40
## integrated_snn_res.2.3 integrated_snn_res.2.4
## AAAGGGCAGTCACGAG--C0 31 32
## AAAGTCCCAGGTTCAT--C0 22 20
## AAAGTCCTCCGAGGCT--C0 22 20
## AACAAAGCACCTGCGA--C0 22 20
## AACCTGACATGCGGTC--C0 31 32
## AACGTCAAGGTTCATC--C0 40 41
## integrated_snn_res.2.5 integrated_snn_res.2.6
## AAAGGGCAGTCACGAG--C0 30 31
## AAAGTCCCAGGTTCAT--C0 23 20
## AAAGTCCTCCGAGGCT--C0 23 20
## AACAAAGCACCTGCGA--C0 23 20
## AACCTGACATGCGGTC--C0 30 31
## AACGTCAAGGTTCATC--C0 41 42
## integrated_snn_res.2.7 integrated_snn_res.2.8
## AAAGGGCAGTCACGAG--C0 33 33
## AAAGTCCCAGGTTCAT--C0 23 22
## AAAGTCCTCCGAGGCT--C0 23 22
## AACAAAGCACCTGCGA--C0 23 22
## AACCTGACATGCGGTC--C0 33 33
## AACGTCAAGGTTCATC--C0 44 44
## integrated_snn_res.2.9 integrated_snn_res.3
## AAAGGGCAGTCACGAG--C0 32 37
## AAAGTCCCAGGTTCAT--C0 22 22
## AAAGTCCTCCGAGGCT--C0 22 22
## AACAAAGCACCTGCGA--C0 22 22
## AACCTGACATGCGGTC--C0 32 37
## AACGTCAAGGTTCATC--C0 43 46
## integrated_snn_res.0.21 integrated_snn_res.0.22
## AAAGGGCAGTCACGAG--C0 6 6
## AAAGTCCCAGGTTCAT--C0 6 6
## AAAGTCCTCCGAGGCT--C0 6 6
## AACAAAGCACCTGCGA--C0 6 6
## AACCTGACATGCGGTC--C0 6 6
## AACGTCAAGGTTCATC--C0 6 6
## integrated_snn_res.0.23 integrated_snn_res.0.24
## AAAGGGCAGTCACGAG--C0 6 7
## AAAGTCCCAGGTTCAT--C0 6 7
## AAAGTCCTCCGAGGCT--C0 6 7
## AACAAAGCACCTGCGA--C0 6 7
## AACCTGACATGCGGTC--C0 6 7
## AACGTCAAGGTTCATC--C0 6 7
## integrated_snn_res.0.25 integrated_snn_res.0.26
## AAAGGGCAGTCACGAG--C0 6 6
## AAAGTCCCAGGTTCAT--C0 6 6
## AAAGTCCTCCGAGGCT--C0 6 6
## AACAAAGCACCTGCGA--C0 6 6
## AACCTGACATGCGGTC--C0 6 6
## AACGTCAAGGTTCATC--C0 6 6
## integrated_snn_res.0.27 integrated_snn_res.0.28
## AAAGGGCAGTCACGAG--C0 6 6
## AAAGTCCCAGGTTCAT--C0 6 6
## AAAGTCCTCCGAGGCT--C0 6 6
## AACAAAGCACCTGCGA--C0 6 6
## AACCTGACATGCGGTC--C0 6 6
## AACGTCAAGGTTCATC--C0 6 6
## integrated_snn_res.0.29 annotation.1 celltype.stim
## AAAGGGCAGTCACGAG--C0 6 Macro Macro_C0
## AAAGTCCCAGGTTCAT--C0 6 Macro Macro_C0
## AAAGTCCTCCGAGGCT--C0 6 Macro Macro_C0
## AACAAAGCACCTGCGA--C0 6 Macro Macro_C0
## AACCTGACATGCGGTC--C0 6 Macro Macro_C0
## AACGTCAAGGTTCATC--C0 6 Macro Macro_C0
## celltype res.0.5.stim annotated.1.sample.id
## AAAGGGCAGTCACGAG--C0 Macro 7_C0 Macro_C0
## AAAGTCCCAGGTTCAT--C0 Macro 7_C0 Macro_C0
## AAAGTCCTCCGAGGCT--C0 Macro 7_C0 Macro_C0
## AACAAAGCACCTGCGA--C0 Macro 7_C0 Macro_C0
## AACCTGACATGCGGTC--C0 Macro 7_C0 Macro_C0
## AACGTCAAGGTTCATC--C0 Macro 7_C0 Macro_C0
## macrophages_subclusters_res.0
## AAAGGGCAGTCACGAG--C0 0
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 0
## AACCTGACATGCGGTC--C0 0
## AACGTCAAGGTTCATC--C0 0
## macrophages_subclusters_res.0.1
## AAAGGGCAGTCACGAG--C0 0
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 0
## AACCTGACATGCGGTC--C0 0
## AACGTCAAGGTTCATC--C0 1
## macrophages_subclusters_res.0.2
## AAAGGGCAGTCACGAG--C0 0
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 0
## AACCTGACATGCGGTC--C0 0
## AACGTCAAGGTTCATC--C0 2
## macrophages_subclusters_res.0.3
## AAAGGGCAGTCACGAG--C0 0
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 0
## AACCTGACATGCGGTC--C0 0
## AACGTCAAGGTTCATC--C0 2
## macrophages_subclusters_res.0.4
## AAAGGGCAGTCACGAG--C0 0
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 0
## AACCTGACATGCGGTC--C0 0
## AACGTCAAGGTTCATC--C0 2
## macrophages_subclusters_res.0.5
## AAAGGGCAGTCACGAG--C0 0
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 0
## AACCTGACATGCGGTC--C0 0
## AACGTCAAGGTTCATC--C0 2
## macrophages_subclusters_res.0.6
## AAAGGGCAGTCACGAG--C0 0
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 0
## AACCTGACATGCGGTC--C0 0
## AACGTCAAGGTTCATC--C0 2
## macrophages_subclusters_res.0.7
## AAAGGGCAGTCACGAG--C0 0
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 0
## AACCTGACATGCGGTC--C0 0
## AACGTCAAGGTTCATC--C0 2
## macrophages_subclusters_res.0.8
## AAAGGGCAGTCACGAG--C0 0
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 0
## AACCTGACATGCGGTC--C0 0
## AACGTCAAGGTTCATC--C0 2
## macrophages_subclusters_res.0.9
## AAAGGGCAGTCACGAG--C0 0
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 0
## AACCTGACATGCGGTC--C0 0
## AACGTCAAGGTTCATC--C0 2
## macrophages_subclusters_res.1
## AAAGGGCAGTCACGAG--C0 0
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 0
## AACCTGACATGCGGTC--C0 0
## AACGTCAAGGTTCATC--C0 2
## macrophages_subclusters_res.1.1
## AAAGGGCAGTCACGAG--C0 0
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 0
## AACCTGACATGCGGTC--C0 0
## AACGTCAAGGTTCATC--C0 2
## macrophages_subclusters_res.1.2
## AAAGGGCAGTCACGAG--C0 3
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 0
## AACCTGACATGCGGTC--C0 0
## AACGTCAAGGTTCATC--C0 2
## macrophages_subclusters_res.1.3
## AAAGGGCAGTCACGAG--C0 2
## AAAGTCCCAGGTTCAT--C0 1
## AAAGTCCTCCGAGGCT--C0 1
## AACAAAGCACCTGCGA--C0 2
## AACCTGACATGCGGTC--C0 2
## AACGTCAAGGTTCATC--C0 3
## macrophages_subclusters_res.1.4
## AAAGGGCAGTCACGAG--C0 2
## AAAGTCCCAGGTTCAT--C0 1
## AAAGTCCTCCGAGGCT--C0 1
## AACAAAGCACCTGCGA--C0 2
## AACCTGACATGCGGTC--C0 2
## AACGTCAAGGTTCATC--C0 3
## macrophages_subclusters_res.1.5
## AAAGGGCAGTCACGAG--C0 2
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 0
## AACAAAGCACCTGCGA--C0 2
## AACCTGACATGCGGTC--C0 2
## AACGTCAAGGTTCATC--C0 3
## macrophages_subclusters_res.1.6
## AAAGGGCAGTCACGAG--C0 3
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 4
## AACAAAGCACCTGCGA--C0 4
## AACCTGACATGCGGTC--C0 3
## AACGTCAAGGTTCATC--C0 2
## macrophages_subclusters_res.1.7
## AAAGGGCAGTCACGAG--C0 3
## AAAGTCCCAGGTTCAT--C0 0
## AAAGTCCTCCGAGGCT--C0 4
## AACAAAGCACCTGCGA--C0 4
## AACCTGACATGCGGTC--C0 3
## AACGTCAAGGTTCATC--C0 1
## macrophages_subclusters_res.1.8
## AAAGGGCAGTCACGAG--C0 3
## AAAGTCCCAGGTTCAT--C0 2
## AAAGTCCTCCGAGGCT--C0 2
## AACAAAGCACCTGCGA--C0 2
## AACCTGACATGCGGTC--C0 3
## AACGTCAAGGTTCATC--C0 7
## macrophages_subclusters_res.1.9
## AAAGGGCAGTCACGAG--C0 3
## AAAGTCCCAGGTTCAT--C0 2
## AAAGTCCTCCGAGGCT--C0 2
## AACAAAGCACCTGCGA--C0 2
## AACCTGACATGCGGTC--C0 8
## AACGTCAAGGTTCATC--C0 7
## macrophages_subclusters_res.2
## AAAGGGCAGTCACGAG--C0 6
## AAAGTCCCAGGTTCAT--C0 2
## AAAGTCCTCCGAGGCT--C0 2
## AACAAAGCACCTGCGA--C0 2
## AACCTGACATGCGGTC--C0 3
## AACGTCAAGGTTCATC--C0 7# install.packages('clustree')
library(clustree)
c <- clustree(macrophages.integrated, prefix = "macrophages_subclusters_res.")
c
pdf(paste0(outdir, "/2_clustree_macrophages.pdf"))
c
dev.off()## png
## 2Using the stability index to assess clusters The stability index from the SC3 package (Kiselev et al. 2017) measures the stability of clusters across resolutions and is automatically calculated when a clustering tree is built.
c <- clustree(macrophages.integrated, prefix = "macrophages_subclusters_res.",
node_colour = "sc3_stability")
pdf(paste0(outdir, "/3_clustree_macrophages_sc3_stability.pdf"))
c
dev.off()## png
## 2SOCS3, IL-1R2, IL1rn
Socs3, Il1r2, Il1rn
c <- clustree(macrophages.integrated, prefix = "macrophages_subclusters_res.",
node_colour = "Socs3", node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/4_clustree_macrophages_Socs3.pdf"))
c
dev.off()## png
## 2c <- clustree(macrophages.integrated, prefix = "macrophages_subclusters_res.",
node_colour = "Il1r2", node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/5_clustree_macrophages_Il1r2.pdf"))
c
dev.off()## png
## 2c <- clustree(macrophages.integrated, prefix = "macrophages_subclusters_res.",
node_colour = "Il1rn", node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/6_clustree_macrophages_Il1rn.pdf"))
c
dev.off()## png
## 2RES 0.3
# install.packages('Polychrome')
library(Polychrome)
set.seed(5658807) # for reproducibility
levels <- levels(as.factor(macrophages.integrated$macrophages_subclusters_res.0.3))
rainbow2.6c <- c("#03539C", "#12999E", "#B7CE05", "#FAEB09",
"#FDA908", "#E82564")
colors.macrophages.subclusters <- createPalette(length(levels),
rainbow2.6c, M = 1000)
names(colors.macrophages.subclusters) <- levels
colors.macrophages.subclusters## 0 1 2 3 <NA> <NA> <NA> <NA>
## "#1662B5" "#0D9A9D" "#C2DB1C" "#FBE22E" "#FCA700" "#F33877" "#FA35F5" "#FFB1DB"
## <NA> <NA>
## "#40F897" "#63691C"slices <- rep(1, length(colors.macrophages.subclusters))
pie(slices, col = colors.macrophages.subclusters, labels = names(colors.macrophages.subclusters))
Idents(macrophages.integrated) <- "macrophages_subclusters_res.0.3"
DimPlot(macrophages.integrated, reduction = "umap", pt.size = 1,
label = T, label.size = 8, cols = colors.macrophages.subclusters)
pdf(paste0(outdir, "/7_DimPlot_umap_macrophages_subclusters_res.0.3.pdf"))
DimPlot(macrophages.integrated, reduction = "umap", pt.size = 0.5,
label = T, label.size = 6, cols = colors.macrophages.subclusters)
dev.off()## png
## 2colors.samples <- c("#12999E", "#FDA908")
names(colors.samples) <- c("C0", "C24")
macrophages.integrated$sample.id <- factor(macrophages.integrated$sample.id,
levels = names(colors.samples))
p <- DimPlot(macrophages.integrated, reduction = "umap", group.by = "sample.id",
cols = colors.samples, label = F)
pdf(paste0(outdir, "/8_UMAP_macrophages_samples.pdf"), width = 10,
height = 9)
p
dev.off()## png
## 2pdf(paste0(outdir, "/9_UMAP_macrophages_samples_noLegend.pdf"),
width = 7, height = 7)
p + NoLegend()
dev.off()## png
## 2Rename subclusters
new.cluster.ids <- paste0("Macro", 1:length(levels(as.factor(macrophages.integrated$macrophages_subclusters_res.0.3))))
names(new.cluster.ids) <- levels(macrophages.integrated)
macrophages.integrated <- RenameIdents(macrophages.integrated,
new.cluster.ids)
macrophages.integrated[["renamed.macrophages.subclusters"]] <- macrophages.integrated@active.identIdents(macrophages.integrated) <- "renamed.macrophages.subclusters"
names(colors.macrophages.subclusters) <- levels(macrophages.integrated$renamed.macrophages.subclusters)
DimPlot(macrophages.integrated, reduction = "umap", pt.size = 1,
label = T, label.size = 8, cols = colors.macrophages.subclusters)
pdf(paste0(outdir, "/10_DimPlot_umap_renamed_macrophages_subclusters.pdf"))
DimPlot(macrophages.integrated, reduction = "umap", pt.size = 0.5,
label = T, label.size = 6, cols = colors.macrophages.subclusters)
dev.off()## png
## 2Idents(macrophages.integrated) <- "renamed.macrophages.subclusters"
DimPlot(macrophages.integrated, reduction = "umap", pt.size = 1,
label = F, label.size = 8, cols = colors.macrophages.subclusters) +
NoLegend()
pdf(paste0(outdir, "/11_DimPlot_umap_renamed_macrophages_subclusters_noLabels_noLegend.pdf"))
DimPlot(macrophages.integrated, reduction = "umap", pt.size = 0.5,
label = F, label.size = 6, cols = colors.macrophages.subclusters) +
NoLegend()
dev.off()## png
## 2saveRDS(macrophages.integrated, paste0(outdir, "/12.macrophages.integrated.rds"))Conserved markers
t <- table(macrophages.integrated$renamed.macrophages.subclusters,
macrophages.integrated$sample.id)
t##
## C0 C24
## Macro1 137 161
## Macro2 72 98
## Macro3 39 50
## Macro4 19 12openxlsx::write.xlsx(as.data.frame.matrix(t), file = paste0(outdir,
"/13_Sample_macrophages_subcluster_composition.xlsx"), rowNames = T,
colNames = T)DefaultAssay(macrophages.integrated) <- "RNA"
Idents(macrophages.integrated) <- "renamed.macrophages.subclusters"
subclusters <- levels(as.factor(macrophages.integrated$renamed.macrophages.subclusters))
conserved.markers <- data.frame(matrix(ncol = 14))
for (c in subclusters) {
print(c)
markers.c <- FindConservedMarkers(macrophages.integrated,
ident.1 = c, grouping.var = "sample.id", verbose = T,
logfc.threshold = -Inf, min.pct = -Inf, min.cells.feature = -Inf,
min.cells.group = -Inf)
markers.c <- cbind(data.frame(cluster = rep(c, dim(markers.c)[1]),
gene = rownames(markers.c)), markers.c)
write.table(markers.c, file = paste0(outdir, "/14_markers_",
c, "_macrophages_subclusters.txt"))
colnames(conserved.markers) <- colnames(markers.c)
conserved.markers <- rbind(conserved.markers, markers.c)
}## [1] "Macro1"
## [1] "Macro2"
## [1] "Macro3"
## [1] "Macro4"# for (c in 1:4) { subcluster <- paste0('AM', c) markers.c <-
# read.table(file = paste0(outdir, '/31_markers_',
# subcluster, '_AM_subclusters.txt'))
# colnames(conserved.markers) <- colnames(markers.c)
# conserved.markers <- rbind(conserved.markers, markers.c) }
# head(conserved.markers)# markers.AM5 <- read.table(file = paste0(outdir,
# '/31_markers_', 'AM5', '_AM_subclusters.txt'))
# colnames(markers.AM5) colnames(conserved.markers)# for ( i in 1:10 ) { markers.AM5[, 34+i] <- rep(NA,
# dim(markers.AM5)[1]) } markers.AM5 <- markers.AM5[, c(1:12,
# 35:44, 13:34)] colnames(markers.AM5) <-
# colnames(conserved.markers) head(markers.AM5)# conserved.markers <- rbind(conserved.markers, markers.AM5)conserved.markers <- conserved.markers[-1, ]
conserved.markers <- conserved.markers[, c(1:2, 13:14, 3:12)]
head(conserved.markers)## cluster gene max_pval minimump_p_val C0_p_val C0_avg_log2FC
## Sdc4 Macro1 Sdc4 6.664112e-02 1.165391e-30 6.664112e-02 0.14916480
## Cd14 Macro1 Cd14 6.988969e-19 5.627031e-30 6.988969e-19 1.71602690
## Cebpb Macro1 Cebpb 1.710274e-03 2.737600e-24 1.710274e-03 0.81625576
## C1qc Macro1 C1qc 4.471371e-17 4.853862e-23 4.471371e-17 1.15050250
## Socs3 Macro1 Socs3 6.343829e-01 3.133599e-22 6.343829e-01 0.09794724
## Pid1 Macro1 Pid1 2.189482e-02 3.136037e-22 2.189482e-02 0.23031958
## C0_pct.1 C0_pct.2 C0_p_val_adj C24_p_val C24_avg_log2FC C24_pct.1
## Sdc4 0.328 0.223 1.000000e+00 5.826955e-31 1.547758 0.994
## Cd14 0.803 0.277 1.565459e-14 2.813516e-30 2.000407 0.981
## Cebpb 0.219 0.092 1.000000e+00 1.368800e-24 1.445991 0.988
## C1qc 0.964 0.546 1.001542e-12 2.426931e-23 1.415284 0.957
## Socs3 0.117 0.100 1.000000e+00 1.566800e-22 1.329294 0.938
## Pid1 0.453 0.331 1.000000e+00 1.568018e-22 1.396260 0.776
## C24_pct.2 C24_p_val_adj
## Sdc4 0.656 1.305180e-26
## Cd14 0.475 6.301994e-26
## Cebpb 0.537 3.065975e-20
## C1qc 0.369 5.436082e-19
## Socs3 0.469 3.509474e-18
## Pid1 0.225 3.512205e-18Only top markers that are positive markers
conserved.markers$marker.type <- apply(conserved.markers, 1, function(x) {
y <- as.numeric(x)
if ( (if (is.na(y[6])) {TRUE} else {y[6]>0})
& (if (is.na(y[11])) {TRUE} else {y[11]>0})
# & (if (is.na(y[16])) {TRUE} else {y[16]>0})
# & (if (is.na(y[21])) {TRUE} else {y[21]>0})
# & (if (is.na(y[26])) {TRUE} else {y[26]>0})
# & (if (is.na(y[31])) {TRUE} else {y[31]>0})
# & (if (is.na(y[36])) {TRUE} else {y[36]>0})
# & (if (is.na(y[41])) {TRUE} else {y[41]>0})
)
{"POS"}
else if ( ( if (is.na(y[6])) {TRUE} else {y[6]<0})
& (if (is.na(y[11])) {TRUE} else {y[11]<0})
# & (if (is.na(y[16])) {TRUE} else {y[16]<0})
# & (if (is.na(y[21])) {TRUE} else {y[21]<0})
# & (if (is.na(y[26])) {TRUE} else {y[26]<0})
# & (if (is.na(y[31])) {TRUE} else {y[31]<0})
# & (if (is.na(y[36])) {TRUE} else {y[36]<0})
# & (if (is.na(y[41])) {TRUE} else {y[41]<0})
)
{"NEG"}
else {"UND"}
})
conserved.markers <- conserved.markers[, c(1:4, 15, 5:14)]openxlsx::write.xlsx(conserved.markers, paste0(outdir, "/15_conserved_markers_macrophages.xlsx"),
colNames = T)
head(conserved.markers)## cluster gene max_pval minimump_p_val marker.type C0_p_val
## Sdc4 Macro1 Sdc4 6.664112e-02 1.165391e-30 POS 6.664112e-02
## Cd14 Macro1 Cd14 6.988969e-19 5.627031e-30 POS 6.988969e-19
## Cebpb Macro1 Cebpb 1.710274e-03 2.737600e-24 POS 1.710274e-03
## C1qc Macro1 C1qc 4.471371e-17 4.853862e-23 POS 4.471371e-17
## Socs3 Macro1 Socs3 6.343829e-01 3.133599e-22 POS 6.343829e-01
## Pid1 Macro1 Pid1 2.189482e-02 3.136037e-22 POS 2.189482e-02
## C0_avg_log2FC C0_pct.1 C0_pct.2 C0_p_val_adj C24_p_val C24_avg_log2FC
## Sdc4 0.14916480 0.328 0.223 1.000000e+00 5.826955e-31 1.547758
## Cd14 1.71602690 0.803 0.277 1.565459e-14 2.813516e-30 2.000407
## Cebpb 0.81625576 0.219 0.092 1.000000e+00 1.368800e-24 1.445991
## C1qc 1.15050250 0.964 0.546 1.001542e-12 2.426931e-23 1.415284
## Socs3 0.09794724 0.117 0.100 1.000000e+00 1.566800e-22 1.329294
## Pid1 0.23031958 0.453 0.331 1.000000e+00 1.568018e-22 1.396260
## C24_pct.1 C24_pct.2 C24_p_val_adj
## Sdc4 0.994 0.656 1.305180e-26
## Cd14 0.981 0.475 6.301994e-26
## Cebpb 0.988 0.537 3.065975e-20
## C1qc 0.957 0.369 5.436082e-19
## Socs3 0.938 0.469 3.509474e-18
## Pid1 0.776 0.225 3.512205e-18saveRDS(macrophages.integrated, paste0(outdir, "/16.macrophages.integrated.rds"))conserved.markers.pos <- conserved.markers[conserved.markers$marker.type ==
"POS", ]
head(conserved.markers.pos)## cluster gene max_pval minimump_p_val marker.type C0_p_val
## Sdc4 Macro1 Sdc4 6.664112e-02 1.165391e-30 POS 6.664112e-02
## Cd14 Macro1 Cd14 6.988969e-19 5.627031e-30 POS 6.988969e-19
## Cebpb Macro1 Cebpb 1.710274e-03 2.737600e-24 POS 1.710274e-03
## C1qc Macro1 C1qc 4.471371e-17 4.853862e-23 POS 4.471371e-17
## Socs3 Macro1 Socs3 6.343829e-01 3.133599e-22 POS 6.343829e-01
## Pid1 Macro1 Pid1 2.189482e-02 3.136037e-22 POS 2.189482e-02
## C0_avg_log2FC C0_pct.1 C0_pct.2 C0_p_val_adj C24_p_val C24_avg_log2FC
## Sdc4 0.14916480 0.328 0.223 1.000000e+00 5.826955e-31 1.547758
## Cd14 1.71602690 0.803 0.277 1.565459e-14 2.813516e-30 2.000407
## Cebpb 0.81625576 0.219 0.092 1.000000e+00 1.368800e-24 1.445991
## C1qc 1.15050250 0.964 0.546 1.001542e-12 2.426931e-23 1.415284
## Socs3 0.09794724 0.117 0.100 1.000000e+00 1.566800e-22 1.329294
## Pid1 0.23031958 0.453 0.331 1.000000e+00 1.568018e-22 1.396260
## C24_pct.1 C24_pct.2 C24_p_val_adj
## Sdc4 0.994 0.656 1.305180e-26
## Cd14 0.981 0.475 6.301994e-26
## Cebpb 0.988 0.537 3.065975e-20
## C1qc 0.957 0.369 5.436082e-19
## Socs3 0.938 0.469 3.509474e-18
## Pid1 0.776 0.225 3.512205e-18table(macrophages.integrated$sample.id)##
## C0 C24
## 267 321colnames(conserved.markers.pos)## [1] "cluster" "gene" "max_pval" "minimump_p_val"
## [5] "marker.type" "C0_p_val" "C0_avg_log2FC" "C0_pct.1"
## [9] "C0_pct.2" "C0_p_val_adj" "C24_p_val" "C24_avg_log2FC"
## [13] "C24_pct.1" "C24_pct.2" "C24_p_val_adj"library(dplyr)
top1.pos <- conserved.markers.pos %>% group_by(cluster) %>% top_n(n = 1,
wt = -minimump_p_val)
top1.pos <- top1.pos %>% group_by(cluster) %>% top_n(n = 1, wt = -max_pval)
top1.pos <- top1.pos %>% group_by(cluster) %>% top_n(n = 1, wt = C24_avg_log2FC)
head(top1.pos)## # A tibble: 4 x 15
## # Groups: cluster [4]
## cluster gene max_pval minimump_p_val marker.type C0_p_val C0_avg_log2FC
## <chr> <chr> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 Macro1 Sdc4 6.66e- 2 1.17e-30 POS 6.66e- 2 0.149
## 2 Macro2 Kap 3.17e-12 1.10e-21 POS 3.17e-12 0.748
## 3 Macro3 Pfkp 2.36e- 2 4.23e-39 POS 2.36e- 2 0.176
## 4 Macro4 Rets… 5.04e- 5 2.83e-21 POS 1.42e-21 0.998
## # … with 8 more variables: C0_pct.1 <dbl>, C0_pct.2 <dbl>, C0_p_val_adj <dbl>,
## # C24_p_val <dbl>, C24_avg_log2FC <dbl>, C24_pct.1 <dbl>, C24_pct.2 <dbl>,
## # C24_p_val_adj <dbl>DefaultAssay(macrophages.integrated) <- "RNA"
f <- FeaturePlot(macrophages.integrated, features = top1.pos$gene,
min.cutoff = "q9", order = T, label = T, label.size = 8,
repel = T)
f
pdf(paste0(outdir, "/17_FeaturePlot_top_pos_markers_macrophages.pdf"),
width = 14, height = 14)
f
dev.off()## png
## 2DotPlot of top markers
DefaultAssay(macrophages.integrated) <- "RNA"
d <- DotPlot(object = macrophages.integrated, features = top1.pos$gene,
cols = c("#03539C", "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.subclusters") +
RotatedAxis()
d
pdf(paste0(outdir, "/18_DotPlot_top_cell_types_pos_markers_clusters.macrophages.pdf"),
width = 7, height = 7)
d
dev.off()## png
## 2DotPlot of top 10 markers
library(dplyr)
top10.pos <- conserved.markers.pos %>% group_by(cluster) %>%
top_n(n = 10, wt = -minimump_p_val)
top10.pos <- top10.pos %>% group_by(cluster) %>% top_n(n = 10,
wt = -max_pval)
top10.pos <- top10.pos %>% group_by(cluster) %>% top_n(n = 10,
wt = C24_avg_log2FC)
head(top10.pos)## # A tibble: 6 x 15
## # Groups: cluster [1]
## cluster gene max_pval minimump_p_val marker.type C0_p_val C0_avg_log2FC
## <chr> <chr> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 Macro1 Sdc4 6.66e- 2 1.17e-30 POS 6.66e- 2 0.149
## 2 Macro1 Cd14 6.99e-19 5.63e-30 POS 6.99e-19 1.72
## 3 Macro1 Cebpb 1.71e- 3 2.74e-24 POS 1.71e- 3 0.816
## 4 Macro1 C1qc 4.47e-17 4.85e-23 POS 4.47e-17 1.15
## 5 Macro1 Socs3 6.34e- 1 3.13e-22 POS 6.34e- 1 0.0979
## 6 Macro1 Pid1 2.19e- 2 3.14e-22 POS 2.19e- 2 0.230
## # … with 8 more variables: C0_pct.1 <dbl>, C0_pct.2 <dbl>, C0_p_val_adj <dbl>,
## # C24_p_val <dbl>, C24_avg_log2FC <dbl>, C24_pct.1 <dbl>, C24_pct.2 <dbl>,
## # C24_p_val_adj <dbl>DefaultAssay(macrophages.integrated) <- "RNA"
d <- DotPlot(object = macrophages.integrated, features = top10.pos$gene,
cols = c("#03539C", "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.subclusters") +
RotatedAxis()
d
pdf(paste0(outdir, "/19_DotPlot_top10_cell_types_pos_markers_clusters.macrophages.pdf"),
width = 14, height = 7)
d
dev.off()## png
## 2Segregation of clusters by various sources of uninteresting variation
Explore additional metrics, such as the number of UMIs and genes per cell, and mitochondrial gene expression by UMAP.
# Determine metrics to plot present in
# tcells.integrated@meta.data
metrics <- c("nCount_RNA", "nFeature_RNA", "percent.mito")
f <- FeaturePlot(macrophages.integrated, reduction = "umap",
features = metrics, pt.size = 0.4, order = TRUE, min.cutoff = "q10",
label = TRUE)
f
pdf(paste0(outdir, "/20_FeaturePlot_macrophages.integrated_umap_metrics.pdf"),
width = 14, height = 14)
f
dev.off()## png
## 2Heatmap
macrophages.integrated <- ScaleData(macrophages.integrated, verbose = T,
features = rownames(macrophages.integrated))
d <- DoHeatmap(macrophages.integrated, features = top10.pos$gene,
group.by = "renamed.macrophages.subclusters", group.colors = colors.macrophages.subclusters) +
NoLegend()
d
pdf(paste0(outdir, "/21_DoHeatmap_top10_pos_macrophages.pdf"),
width = 15, height = 13)
d
dev.off()## png
## 2DotPlot with known/useful markers
immune.cells <- c("Ptprc", "Spi1", "Cd3g")
b.cells <- c("Cd19", "Cd79a", "Cd79b")
t.cells <- "Cd3e"
nk.cells <- c("Ncr1")
myeloid.cells <- "Cd68"
neutrophils <- c("Csf3r", "S100a8", "S100a9", "Retnlg")
monocytes.macrophages.dendritic.cells <- "Csf1r"
dc <- c("Siglech", "Cd209a")
dc1 <- "Xcr1"
dc2 <- c("Ccl18", "Nr4a3", "Tnfsf12")
dc3 <- "Fscn1"
pdc <- "Tcf4"
monocytes <- c("Ly6c1", "Ly6c2", "Ccr2", "Cx3cr1", "Fcgr4")
macrophages <- c("Mrc1", "Lyz2", "Adgre1", "Axl", "Cxcl2", "Trem2",
"C1qc", "Fcgr1", "C5ar1")
am.markers <- c("Pparg", "Itgax", "Ear1", "Ear2", "Ear3", "Car4",
"Siglec5", "SiglecF")
im.markers <- c("Itgam", "Cx3cr1", "Ccl12", "Bhlme40", "Trem2",
"Cxcl12", "Csf1")
peribronchial.macrophages <- c("Ccr2", "Bhlhe40", "Bhlhe41")
perivascular.macrophages <- c("Lyve1", "F13a1")
proliferating.cells <- "Mki67"
mesenchymal.cells <- c("Col1a1", "Col1a2", "Pdgfra")
epithelial.cells <- c("Epcam", "Cdh1")
mast.cells <- "Mcpt8"
endothelial.cells <- c("Pecam1", "Emcn")markers <- unique(c(immune.cells, b.cells, t.cells, nk.cells,
myeloid.cells, neutrophils, monocytes.macrophages.dendritic.cells,
dc, dc1, dc2, dc3, pdc, monocytes, macrophages, am.markers,
im.markers, peribronchial.macrophages, perivascular.macrophages,
proliferating.cells, mesenchymal.cells, epithelial.cells,
mast.cells, endothelial.cells))Dot plot
DefaultAssay(macrophages.integrated) <- "RNA"
d <- DotPlot(object = macrophages.integrated, features = markers,
cols = c("#03539C", "#E82564"), dot.scale = 15, group.by = "renamed.macrophages.subclusters") +
RotatedAxis()
d
pdf(paste0(outdir, "/22_DotPlot_known_markers_macrophages.pdf"),
width = 18, height = 4)
d
dev.off()## png
## 2Expression of genes of interest: Socs3, Il1r2, Il1rn
goi <- c("Socs3", "Il1r2", "Il1rn")DefaultAssay(macrophages.integrated) <- "RNA"
f <- FeaturePlot(macrophages.integrated, features = goi, min.cutoff = "q9",
order = T, label = T, label.size = 8, repel = T)
f
pdf(paste0(outdir, "/23_FeaturePlot_goi_macrophages.pdf"), width = 14,
height = 14)
f
dev.off()## png
## 2DefaultAssay(macrophages.integrated) <- "RNA"
d <- DotPlot(object = macrophages.integrated, features = goi,
cols = c("#03539C", "#E82564"), dot.scale = 15, group.by = "renamed.macrophages.subclusters") +
RotatedAxis()
d
pdf(paste0(outdir, "/24_DotPlot_goi_macrophages.pdf"), width = 6,
height = 4)
d
dev.off()## png
## 2DIFFERENTIAL EXPRESSION ANALYSES
differential expression studies for the following groups: 1. C24 vs C0 within each cell type (from manual annotation annotation.1)
1. C24 vs C0 within each cell type (from manual annotation annotation.1)
macrophages.integrated$renamed.macrophages.subclusters.sample.id <- paste(macrophages.integrated$renamed.macrophages.subclusters,
macrophages.integrated$sample.id, sep = "_")
levels(factor(macrophages.integrated$renamed.macrophages.subclusters.sample.id))## [1] "Macro1_C0" "Macro1_C24" "Macro2_C0" "Macro2_C24" "Macro3_C0"
## [6] "Macro3_C24" "Macro4_C0" "Macro4_C24"DefaultAssay(macrophages.integrated) <- "RNA"
output <- paste0(outdir, "/25_DE_macrophages_subclusters_C24_vs_C0_")
levels(factor(macrophages.integrated$renamed.macrophages.subclusters.sample.id))## [1] "Macro1_C0" "Macro1_C24" "Macro2_C0" "Macro2_C24" "Macro3_C0"
## [6] "Macro3_C24" "Macro4_C0" "Macro4_C24"for (cell.type in levels(factor(macrophages.integrated$renamed.macrophages.subclusters))) {
try({
print(cell.type)
ident1 <- paste0(cell.type, "_C24")
ident2 <- paste0(cell.type, "_C0")
Idents(macrophages.integrated) <- "renamed.macrophages.subclusters.sample.id"
condition.diffgenes <- FindMarkers(macrophages.integrated,
ident.1 = ident1, ident.2 = ident2, logfc.threshold = -Inf,
test.use = "wilcox", min.pct = -Inf, verbose = T,
min.cells.feature = -Inf, min.cells.group = -Inf)
condition.diffgenes.adjpval.0.05 <- condition.diffgenes[condition.diffgenes$p_val_adj <
0.05, ]
condition.diffgenes.adjpval.0.05.upregulated.sorted <- condition.diffgenes.adjpval.0.05[condition.diffgenes.adjpval.0.05$avg_log2FC >
0, ]
condition.diffgenes.adjpval.0.05.upregulated.sorted <- condition.diffgenes.adjpval.0.05.upregulated.sorted[order(condition.diffgenes.adjpval.0.05.upregulated.sorted$avg_log2FC,
decreasing = T), ]
condition.diffgenes.adjpval.0.05.downregulated.sorted <- condition.diffgenes.adjpval.0.05[condition.diffgenes.adjpval.0.05$avg_log2FC <
0, ]
condition.diffgenes.adjpval.0.05.downregulated.sorted <- condition.diffgenes.adjpval.0.05.downregulated.sorted[order(condition.diffgenes.adjpval.0.05.downregulated.sorted$avg_log2FC,
decreasing = F), ]
list_of_datasets <- list(sig.upregulated.genes.ordered.by.avg.logFC = condition.diffgenes.adjpval.0.05.upregulated.sorted,
sig.downregulated.genes.ordered.by.avg.logFC = condition.diffgenes.adjpval.0.05.downregulated.sorted,
all.genes.ordered.by.adj.pval = condition.diffgenes)
openxlsx::write.xlsx(list_of_datasets, file = paste0(output,
cell.type, ".xlsx"), row.names = T)
})
}## [1] "Macro1"
## [1] "Macro2"
## [1] "Macro3"
## [1] "Macro4"for (cell.type in levels(factor(macrophages.integrated$renamed.macrophages.subclusters))) {
try({
print(cell.type)
filename <- paste0(output, cell.type, ".xlsx")
sheets <- openxlsx::getSheetNames(filename)
SheetList <- lapply(sheets, openxlsx::read.xlsx, xlsxFile = filename)
names(SheetList) <- sheets
condition.diffgenes.adjpval.0.05.upregulated.sorted <- SheetList[[1]]
rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted) <- condition.diffgenes.adjpval.0.05.upregulated.sorted[,
1]
condition.diffgenes.adjpval.0.05.upregulated.sorted <- condition.diffgenes.adjpval.0.05.upregulated.sorted[,
-1]
condition.diffgenes.adjpval.0.05.downregulated.sorted <- SheetList[[2]]
rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted) <- condition.diffgenes.adjpval.0.05.downregulated.sorted[,
1]
condition.diffgenes.adjpval.0.05.downregulated.sorted <- condition.diffgenes.adjpval.0.05.downregulated.sorted[,
-1]
Idents(annotated) <- "renamed.macrophages.subclusters"
try({
f <- FeaturePlot(macrophages.integrated, features = rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted)[1],
split.by = "sample.id", cols = c("grey", "#E82564"),
order = T, min.cutoff = "q10", max.cutoff = "q90")
pdf(paste0(output, cell.type, "_FeaturePlot_top_upregulated_split_by_sample.pdf"),
width = 14, height = 7)
print(f)
dev.off()
})
try({
f <- FeaturePlot(macrophages.integrated, features = rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted)[1],
split.by = "sample.id", cols = c("grey", "#E82564"),
order = T, min.cutoff = "q10", max.cutoff = "q90")
pdf(paste0(output, cell.type, "_FeaturePlot_top_downregulated_split_by_sample.pdf"),
width = 14, height = 7)
print(f)
dev.off()
})
try({
plot <- VlnPlot(macrophages.integrated, features = rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted)[1],
group.by = "renamed.macrophages.subclusters",
split.by = "sample.id", pt.size = 0, combine = FALSE,
cols = colors.samples)
pdf(paste0(output, cell.type, "_VlnPlot_top_upregulated_group_by_cluster_split_by_sample.pdf"),
width = 11, height = 8.5)
print(plot)
dev.off()
})
try({
plot <- VlnPlot(macrophages.integrated, features = rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted)[1],
group.by = "renamed.macrophages.subclusters",
split.by = "sample.id", pt.size = 0, combine = FALSE,
cols = colors.samples)
pdf(paste0(output, cell.type, "_VlnPlot_top_downregulated_group_by_cluster_split_by_sample.pdf"),
width = 11, height = 8.5)
print(plot)
dev.off()
})
## DotPlot with 50 most upregulated genes
Idents(macrophages.integrated) <- "renamed.macrophages.subclusters.sample.id"
try({
d <- DotPlot(macrophages.integrated, features = rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted)[1:50],
cols = c("#03539C", "#E82564"), dot.scale = 8,
group.by = "renamed.macrophages.subclusters.sample.id") +
RotatedAxis()
pdf(paste0(output, cell.type, "_DotPlot_top50_upregulated_groupby_annotation.1.sample.id.pdf"),
width = (4 + min(50, length(rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted))) *
15/50), height = 17)
print(d)
dev.off()
})
## DotPlot with 50 most downregulated genes
Idents(macrophages.integrated) <- "renamed.macrophages.subclusters.sample.id"
try({
d <- DotPlot(macrophages.integrated, features = rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted)[1:50],
cols = c("#03539C", "#E82564"), dot.scale = 8,
group.by = "renamed.macrophages.subclusters.sample.id") +
RotatedAxis()
pdf(paste0(output, cell.type, "_DotPlot_top50_downregulated_groupby_annotation.1.sample.id.pdf"),
width = (4 + min(50, length(rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted))) *
15/50), height = 17)
print(d)
dev.off()
})
})
}## [1] "Macro1"## [[1]]##
## [[1]]## [1] "Macro2"## [[1]]##
## [[1]]## [1] "Macro3"## [[1]]##
## [[1]]## [1] "Macro4"## Error : None of the requested features were found: NA in slot data
## [[1]]##
## Error in FetchData(object = object, vars = features, slot = slot) :
## None of the requested variables were found: NA## Error in h(simpleError(msg, call)) :
## error in evaluating the argument 'x' in selecting a method for function 'mean': non-numeric argument to mathematical functionPlots for genes of interest
for (gene in goi) {
try({
Idents(macrophages.integrated) <- "renamed.macrophages.subclusters"
try({
f <- FeaturePlot(macrophages.integrated, features = gene,
split.by = "sample.id", cols = c("grey", "#E82564"),
order = T, min.cutoff = "q10", max.cutoff = "q90",
label = T)
pdf(paste0(outdir, "/26_", gene, "_FeaturePlot_split_by_sample.pdf"),
width = 14, height = 7)
print(f)
dev.off()
})
try({
plot <- VlnPlot(macrophages.integrated, features = gene,
group.by = "renamed.macrophages.subclusters",
split.by = "sample.id", pt.size = 0, combine = FALSE,
cols = colors.samples)
pdf(paste0(outdir, "/27_", gene, "_VlnPlot_group_by_cluster_split_by_sample.pdf"),
width = 11, height = 8.5)
print(plot)
dev.off()
})
})
}## [[1]]## [[1]]## [[1]]Idents(macrophages.integrated) <- "renamed.macrophages.subclusters.sample.id"
try({
d <- DotPlot(macrophages.integrated, features = goi, cols = c("#03539C",
"#E82564"), dot.scale = 8, group.by = "renamed.macrophages.subclusters.sample.id") +
RotatedAxis()
pdf(paste0(outdir, "/28_", paste(goi, sep = "_"), "_DotPlot_groupby_annotation.1.sample.id.pdf"),
width = (4 + length(goi) * 15/50), height = 17)
print(d)
dev.off()
})## png
## 2Plots for more genes of interest: F480 (antibody), CD11b, CD11c, CD64, MHCII
In mouse: F4/80 (antibody) > Emr1, Adgre1 CD11b -> Itgam CD11c -> Itgax CD64 -> Fcgr1 MHCII -> H2-Ab
goi <- c("Adgre1", "Itgam", "Itgax", "Fcgr1", "H2-Ab1")DefaultAssay(macrophages.integrated) <- "RNA"
for (gene in goi) {
c <- clustree(macrophages.integrated, prefix = "macrophages_subclusters_res.",
node_colour = gene, node_colour_aggr = "mean")
pdf(paste0(outdir, "/29_clustree_", gene, "_macrophages.pdf"))
print(c)
dev.off()
}for (gene in goi) {
try({
Idents(macrophages.integrated) <- "renamed.macrophages.subclusters"
try({
f <- FeaturePlot(macrophages.integrated, features = gene,
split.by = "sample.id", cols = c("grey", "#E82564"),
order = T, min.cutoff = "q10", max.cutoff = "q90",
label = T)
pdf(paste0(outdir, "/30_", gene, "_FeaturePlot_split_by_sample.pdf"),
width = 14, height = 7)
print(f)
dev.off()
})
try({
plot <- VlnPlot(macrophages.integrated, features = gene,
group.by = "renamed.macrophages.subclusters",
split.by = "sample.id", pt.size = 0, combine = FALSE,
cols = colors.samples)
pdf(paste0(outdir, "/31_", gene, "_VlnPlot_group_by_cluster_split_by_sample.pdf"),
width = 11, height = 8.5)
print(plot)
dev.off()
})
})
}## [[1]]## [[1]]## [[1]]## [[1]]## [[1]]Idents(macrophages.integrated) <- "renamed.macrophages.subclusters.sample.id"
try({
d <- DotPlot(macrophages.integrated, features = goi, cols = c("#03539C",
"#E82564"), dot.scale = 8, group.by = "renamed.macrophages.subclusters.sample.id") +
RotatedAxis()
pdf(paste0(outdir, "/32_", paste(goi, sep = "_"), "_DotPlot_groupby_macrophages_subclusters.sample.id.pdf"),
width = (4 + length(goi) * 15/50), height = 17)
print(d)
dev.off()
})## png
## 2Expression of marker Tim3/Havcr2 from reference Hasegawa 2021
https://jasn.asnjournals.org/content/early/2021/04/17/ASN.2020121723
DefaultAssay(macrophages.integrated) <- "RNA"
Idents(macrophages.integrated) <- "renamed.macrophages.subclusters"
f <- FeaturePlot(macrophages.integrated, features = "Havcr2",
min.cutoff = "q9", order = T, label = T, label.size = 8,
repel = T)
f
pdf(paste0(outdir, "/33_FeaturePlot_Havcr2_macrophages.pdf"),
width = 7, height = 7)
f
dev.off()## png
## 2QC metrics to possibly remove high mitochondrial expression cells
v <- VlnPlot(macrophages.integrated, features = c("nFeature_RNA",
"nCount_RNA", "percent.mito"), ncol = 3, cols = colors.samples,
group.by = "sample.id")
v
pdf(paste0(outdir, "/34_Seurat_QC_VlnPlot_samples.pdf"), width = 21,
height = 7)
v
dev.off()## png
## 2summary(macrophages.integrated$percent.mito)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.7848 4.7934 7.9638 12.3109 16.1364 54.687820% seems like a reasonable threshold Jamie: Regarding your other question, I do think having a lower threshold for mitochondrial genes would be good in the macrophage subpopulations. I think the high mitochondrial content is mainly applicable to the tubular epithelial clusters.
First clean up of the subset by removing:
- genes that are not expressed in any cell
dim(macrophages.integrated)## [1] 22399 588keep_feature <- rowSums(as.matrix(GetAssayData(macrophages.integrated,
slot = "counts"))) > 0
summary(keep_feature)## Mode FALSE TRUE
## logical 8106 14293- cells in which no gene expression is detected
keep_cell <- colSums(as.matrix(GetAssayData(macrophages.integrated,
slot = "counts"))) > 0
summary(keep_cell)## Mode TRUE
## logical 588macrophages.integrated.filtered <- macrophages.integrated[keep_feature,
keep_cell]
dim(macrophages.integrated.filtered)## [1] 14293 588Remove cells with mito expn > 20%
dim(macrophages.integrated.filtered)## [1] 14293 588keep_cell <- macrophages.integrated.filtered$percent.mito <=
20
summary(keep_cell)## Mode FALSE TRUE
## logical 108 480macrophages.integrated.filtered <- macrophages.integrated.filtered[,
keep_cell]
dim(macrophages.integrated.filtered)## [1] 14293 480New subclsutering after filtering
DefaultAssay(macrophages.integrated.filtered) <- "RNA"
macrophages.2.list <- SplitObject(macrophages.integrated.filtered,
split.by = "sample.id")
macrophages.2.list## $C0
## An object of class Seurat
## 16293 features across 212 samples within 2 assays
## Active assay: RNA (14293 features, 0 variable features)
## 1 other assay present: integrated
## 2 dimensional reductions calculated: pca, umap
##
## $C24
## An object of class Seurat
## 16293 features across 268 samples within 2 assays
## Active assay: RNA (14293 features, 0 variable features)
## 1 other assay present: integrated
## 2 dimensional reductions calculated: pca, umapfor (i in 1:length(macrophages.2.list)) {
macrophages.2.list[[i]] <- FindVariableFeatures(macrophages.2.list[[i]],
selection.method = "vst", nfeatures = 2000, verbose = FALSE)
}k.filter <- min(200, min(sapply(macrophages.2.list, ncol)))
k.filter## [1] 200macrophages.2.anchors <- FindIntegrationAnchors(object.list = macrophages.2.list,
dims = 1:30, k.filter = k.filter)macrophages.2.integrated <- IntegrateData(anchorset = macrophages.2.anchors,
dims = 1:30)DefaultAssay(macrophages.2.integrated) <- "integrated"
macrophages.2.integrated <- ScaleData(macrophages.2.integrated,
verbose = T)dim(macrophages.2.integrated)## [1] 2000 480macrophages.2.integrated <- RunPCA(macrophages.2.integrated,
npcs = 50, verbose = T)system.time(macrophages.2.integrated <- JackStraw(macrophages.2.integrated,
num.replicate = 100, dims = 50))## user system elapsed
## 63.442 0.093 63.595macrophages.2.integrated <- ScoreJackStraw(macrophages.2.integrated,
dims = 1:50)
j <- JackStrawPlot(macrophages.2.integrated, dims = 1:50)
pdf(paste0(outdir, "/35_JackStrawPlot_macrophages_2.pdf"))
j
dev.off()## png
## 2macrophages.2.integrated <- FindNeighbors(macrophages.2.integrated,
dims = 1:48)
macrophages.2.integrated <- RunUMAP(macrophages.2.integrated,
dims = 1:48)for (i in seq(0, 2, 0.1)) {
macrophages.2.integrated <- FindClusters(macrophages.2.integrated,
resolution = i, verbose = F)
macrophages.2.integrated[[paste0("macrophages_2_subclusters_res.",
i)]] <- macrophages.2.integrated$seurat_clusters
}# install.packages('clustree')
library(clustree)
c <- clustree(macrophages.2.integrated, prefix = "macrophages_2_subclusters_res.")
c
pdf(paste0(outdir, "/36_clustree_macrophages_2.pdf"))
c
dev.off()## png
## 2Using the stability index to assess clusters The stability index from the SC3 package (Kiselev et al. 2017) measures the stability of clusters across resolutions and is automatically calculated when a clustering tree is built.
c <- clustree(macrophages.2.integrated, prefix = "macrophages_2_subclusters_res.",
node_colour = "sc3_stability")
pdf(paste0(outdir, "/37_clustree_macrophages_2_sc3_stability.pdf"))
c
dev.off()## png
## 2SOCS3, IL-1R2, IL1rn
Socs3, Il1r2, Il1rn
c <- clustree(macrophages.2.integrated, prefix = "macrophages_2_subclusters_res.",
node_colour = "Socs3", node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/38_clustree_macrophages_2_Socs3.pdf"))
c
dev.off()## png
## 2c <- clustree(macrophages.2.integrated, prefix = "macrophages_2_subclusters_res.",
node_colour = "Il1r2", node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/39_clustree_macrophages_2_Il1r2.pdf"))
c
dev.off()## png
## 2c <- clustree(macrophages.2.integrated, prefix = "macrophages_2_subclusters_res.",
node_colour = "Il1rn", node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/40_clustree_macrophages_2_Il1rn.pdf"))
c
dev.off()## png
## 2RES 0.2
# install.packages('Polychrome')
library(Polychrome)
set.seed(5658807) # for reproducibility
levels <- levels(as.factor(macrophages.2.integrated$macrophages_2_subclusters_res.0.2))
rainbow2.6c <- c("#03539C", "#12999E", "#B7CE05", "#FAEB09",
"#FDA908", "#E82564")
colors.macrophages.2.subclusters <- createPalette(length(levels),
rainbow2.6c, M = 1000)
names(colors.macrophages.2.subclusters) <- levels
colors.macrophages.2.subclusters## 0 1 2 <NA> <NA> <NA> <NA> <NA>
## "#1662B5" "#0D9A9D" "#C2DB1C" "#FBE22E" "#FCA700" "#F33877" "#FA35F5" "#FFB1DB"
## <NA> <NA> <NA>
## "#40F897" "#63691C" "#4F00FF"slices <- rep(1, length(colors.macrophages.2.subclusters))
pie(slices, col = colors.macrophages.2.subclusters, labels = names(colors.macrophages.2.subclusters))
Idents(macrophages.2.integrated) <- "macrophages_2_subclusters_res.0.2"
DimPlot(macrophages.2.integrated, reduction = "umap", pt.size = 1,
label = T, label.size = 8, cols = colors.macrophages.2.subclusters)
pdf(paste0(outdir, "/41_DimPlot_umap_macrophages_2_subclusters_res.0.2.pdf"))
DimPlot(macrophages.2.integrated, reduction = "umap", pt.size = 1,
label = T, label.size = 6, cols = colors.macrophages.2.subclusters)
dev.off()## png
## 2colors.samples <- c("#12999E", "#FDA908")
names(colors.samples) <- c("C0", "C24")
macrophages.2.integrated$sample.id <- factor(macrophages.2.integrated$sample.id,
levels = names(colors.samples))
p <- DimPlot(macrophages.2.integrated, reduction = "umap", group.by = "sample.id",
cols = colors.samples, label = F)
pdf(paste0(outdir, "/42_UMAP_macrophages_2_samples.pdf"), width = 10,
height = 9)
p
dev.off()## png
## 2pdf(paste0(outdir, "/43_UMAP_macrophages_2_samples_noLegend.pdf"),
width = 7, height = 7)
p + NoLegend()
dev.off()## png
## 2Rename subclusters
new.cluster.ids <- paste0("Macro", 1:length(levels(as.factor(macrophages.2.integrated$macrophages_2_subclusters_res.0.2))))
names(new.cluster.ids) <- levels(macrophages.2.integrated)
macrophages.2.integrated <- RenameIdents(macrophages.2.integrated,
new.cluster.ids)
macrophages.2.integrated[["renamed.macrophages.2.subclusters"]] <- macrophages.2.integrated@active.identIdents(macrophages.2.integrated) <- "renamed.macrophages.2.subclusters"
names(colors.macrophages.2.subclusters) <- levels(macrophages.2.integrated$renamed.macrophages.2.subclusters)
DimPlot(macrophages.2.integrated, reduction = "umap", pt.size = 1,
label = T, label.size = 8, cols = colors.macrophages.2.subclusters)
pdf(paste0(outdir, "/44_DimPlot_umap_renamed_macrophages_2_subclusters.pdf"))
DimPlot(macrophages.2.integrated, reduction = "umap", pt.size = 1,
label = T, label.size = 6, cols = colors.macrophages.2.subclusters)
dev.off()## png
## 2Idents(macrophages.2.integrated) <- "renamed.macrophages.2.subclusters"
DimPlot(macrophages.2.integrated, reduction = "umap", pt.size = 1,
label = F, label.size = 8, cols = colors.macrophages.2.subclusters) +
NoLegend()
pdf(paste0(outdir, "/45_DimPlot_umap_renamed_macrophages_2_subclusters_noLabels_noLegend.pdf"))
DimPlot(macrophages.2.integrated, reduction = "umap", pt.size = 1,
label = F, label.size = 6, cols = colors.macrophages.2.subclusters) +
NoLegend()
dev.off()## png
## 2saveRDS(macrophages.2.integrated, paste0(outdir, "/46.macrophages.2.integrated.rds"))Conserved markers
t <- table(macrophages.2.integrated$renamed.macrophages.2.subclusters,
macrophages.2.integrated$sample.id)
t##
## C0 C24
## Macro1 139 151
## Macro2 38 75
## Macro3 35 42openxlsx::write.xlsx(as.data.frame.matrix(t), file = paste0(outdir,
"/47_Sample_macrophages_2_subcluster_composition.xlsx"),
rowNames = T, colNames = T)DefaultAssay(macrophages.2.integrated) <- "RNA"
Idents(macrophages.2.integrated) <- "renamed.macrophages.2.subclusters"
subclusters <- levels(as.factor(macrophages.2.integrated$renamed.macrophages.2.subclusters))
conserved.markers <- data.frame(matrix(ncol = 14))
for (c in subclusters) {
print(c)
markers.c <- FindConservedMarkers(macrophages.2.integrated,
ident.1 = c, grouping.var = "sample.id", verbose = T,
logfc.threshold = -Inf, min.pct = -Inf, min.cells.feature = -Inf,
min.cells.group = -Inf)
markers.c <- cbind(data.frame(cluster = rep(c, dim(markers.c)[1]),
gene = rownames(markers.c)), markers.c)
write.table(markers.c, file = paste0(outdir, "/48_markers_",
c, "_macrophages_2_subclusters.txt"))
colnames(conserved.markers) <- colnames(markers.c)
conserved.markers <- rbind(conserved.markers, markers.c)
}## [1] "Macro1"
## [1] "Macro2"
## [1] "Macro3"# for (c in 1:4) { subcluster <- paste0('AM', c) markers.c <-
# read.table(file = paste0(outdir, '/31_markers_',
# subcluster, '_AM_subclusters.txt'))
# colnames(conserved.markers) <- colnames(markers.c)
# conserved.markers <- rbind(conserved.markers, markers.c) }
# head(conserved.markers)# markers.AM5 <- read.table(file = paste0(outdir,
# '/31_markers_', 'AM5', '_AM_subclusters.txt'))
# colnames(markers.AM5) colnames(conserved.markers)# for ( i in 1:10 ) { markers.AM5[, 34+i] <- rep(NA,
# dim(markers.AM5)[1]) } markers.AM5 <- markers.AM5[, c(1:12,
# 35:44, 13:34)] colnames(markers.AM5) <-
# colnames(conserved.markers) head(markers.AM5)# conserved.markers <- rbind(conserved.markers, markers.AM5)conserved.markers <- conserved.markers[-1, ]
conserved.markers <- conserved.markers[, c(1:2, 13:14, 3:12)]
head(conserved.markers)## cluster gene max_pval minimump_p_val C0_p_val C0_avg_log2FC
## Sdc4 Macro1 Sdc4 2.202498e-01 9.940550e-25 2.202498e-01 -0.07585808
## Cd14 Macro1 Cd14 3.933950e-13 2.413485e-23 3.933950e-13 1.97725880
## Cebpb Macro1 Cebpb 1.065972e-02 1.624807e-20 1.065972e-02 0.66786332
## Pla2g7 Macro1 Pla2g7 1.069071e-05 1.458113e-18 1.069071e-05 1.31291835
## Pid1 Macro1 Pid1 1.655140e-01 4.899195e-18 1.655140e-01 0.04600540
## Ctsb Macro1 Ctsb 1.454003e-06 1.396807e-17 1.454003e-06 0.84584724
## C0_pct.1 C0_pct.2 C0_p_val_adj C24_p_val C24_avg_log2FC C24_pct.1
## Sdc4 0.309 0.219 1.000000e+00 4.970275e-25 1.494107 1.000
## Cd14 0.755 0.233 5.622794e-09 1.206743e-23 1.878361 0.980
## Cebpb 0.194 0.068 1.000000e+00 8.124037e-21 1.492180 0.987
## Pla2g7 0.324 0.055 1.528024e-01 7.290567e-19 1.376837 0.709
## Pid1 0.460 0.384 1.000000e+00 2.449597e-18 1.276004 0.788
## Ctsb 0.806 0.534 2.078207e-02 6.984035e-18 1.269536 0.967
## C24_pct.2 C24_p_val_adj
## Sdc4 0.667 7.104014e-21
## Cd14 0.496 1.724797e-19
## Cebpb 0.556 1.161169e-16
## Pla2g7 0.103 1.042041e-14
## Pid1 0.231 3.501209e-14
## Ctsb 0.556 9.982281e-14Only top markers that are positive markers
conserved.markers$marker.type <- apply(conserved.markers, 1, function(x) {
y <- as.numeric(x)
if ( (if (is.na(y[6])) {TRUE} else {y[6]>0})
& (if (is.na(y[11])) {TRUE} else {y[11]>0})
# & (if (is.na(y[16])) {TRUE} else {y[16]>0})
# & (if (is.na(y[21])) {TRUE} else {y[21]>0})
# & (if (is.na(y[26])) {TRUE} else {y[26]>0})
# & (if (is.na(y[31])) {TRUE} else {y[31]>0})
# & (if (is.na(y[36])) {TRUE} else {y[36]>0})
# & (if (is.na(y[41])) {TRUE} else {y[41]>0})
)
{"POS"}
else if ( ( if (is.na(y[6])) {TRUE} else {y[6]<0})
& (if (is.na(y[11])) {TRUE} else {y[11]<0})
# & (if (is.na(y[16])) {TRUE} else {y[16]<0})
# & (if (is.na(y[21])) {TRUE} else {y[21]<0})
# & (if (is.na(y[26])) {TRUE} else {y[26]<0})
# & (if (is.na(y[31])) {TRUE} else {y[31]<0})
# & (if (is.na(y[36])) {TRUE} else {y[36]<0})
# & (if (is.na(y[41])) {TRUE} else {y[41]<0})
)
{"NEG"}
else {"UND"}
})
conserved.markers <- conserved.markers[, c(1:4, 15, 5:14)]openxlsx::write.xlsx(conserved.markers, paste0(outdir, "/49_conserved_markers_macrophages_2.xlsx"),
colNames = T)
head(conserved.markers)## cluster gene max_pval minimump_p_val marker.type C0_p_val
## Sdc4 Macro1 Sdc4 2.202498e-01 9.940550e-25 UND 2.202498e-01
## Cd14 Macro1 Cd14 3.933950e-13 2.413485e-23 POS 3.933950e-13
## Cebpb Macro1 Cebpb 1.065972e-02 1.624807e-20 POS 1.065972e-02
## Pla2g7 Macro1 Pla2g7 1.069071e-05 1.458113e-18 POS 1.069071e-05
## Pid1 Macro1 Pid1 1.655140e-01 4.899195e-18 POS 1.655140e-01
## Ctsb Macro1 Ctsb 1.454003e-06 1.396807e-17 POS 1.454003e-06
## C0_avg_log2FC C0_pct.1 C0_pct.2 C0_p_val_adj C24_p_val C24_avg_log2FC
## Sdc4 -0.07585808 0.309 0.219 1.000000e+00 4.970275e-25 1.494107
## Cd14 1.97725880 0.755 0.233 5.622794e-09 1.206743e-23 1.878361
## Cebpb 0.66786332 0.194 0.068 1.000000e+00 8.124037e-21 1.492180
## Pla2g7 1.31291835 0.324 0.055 1.528024e-01 7.290567e-19 1.376837
## Pid1 0.04600540 0.460 0.384 1.000000e+00 2.449597e-18 1.276004
## Ctsb 0.84584724 0.806 0.534 2.078207e-02 6.984035e-18 1.269536
## C24_pct.1 C24_pct.2 C24_p_val_adj
## Sdc4 1.000 0.667 7.104014e-21
## Cd14 0.980 0.496 1.724797e-19
## Cebpb 0.987 0.556 1.161169e-16
## Pla2g7 0.709 0.103 1.042041e-14
## Pid1 0.788 0.231 3.501209e-14
## Ctsb 0.967 0.556 9.982281e-14saveRDS(macrophages.2.integrated, paste0(outdir, "/50.macrophages.2.integrated.rds"))conserved.markers.pos <- conserved.markers[conserved.markers$marker.type ==
"POS", ]
head(conserved.markers.pos)## cluster gene max_pval minimump_p_val marker.type C0_p_val
## Cd14 Macro1 Cd14 3.933950e-13 2.413485e-23 POS 3.933950e-13
## Cebpb Macro1 Cebpb 1.065972e-02 1.624807e-20 POS 1.065972e-02
## Pla2g7 Macro1 Pla2g7 1.069071e-05 1.458113e-18 POS 1.069071e-05
## Pid1 Macro1 Pid1 1.655140e-01 4.899195e-18 POS 1.655140e-01
## Ctsb Macro1 Ctsb 1.454003e-06 1.396807e-17 POS 1.454003e-06
## C1qa Macro1 C1qa 5.312464e-13 1.480332e-16 POS 7.401660e-17
## C0_avg_log2FC C0_pct.1 C0_pct.2 C0_p_val_adj C24_p_val C24_avg_log2FC
## Cd14 1.9772588 0.755 0.233 5.622794e-09 1.206743e-23 1.878361
## Cebpb 0.6678633 0.194 0.068 1.000000e+00 8.124037e-21 1.492180
## Pla2g7 1.3129183 0.324 0.055 1.528024e-01 7.290567e-19 1.376837
## Pid1 0.0460054 0.460 0.384 1.000000e+00 2.449597e-18 1.276004
## Ctsb 0.8458472 0.806 0.534 2.078207e-02 6.984035e-18 1.269536
## C1qa 1.6955557 0.957 0.452 1.057919e-12 5.312464e-13 1.032159
## C24_pct.1 C24_pct.2 C24_p_val_adj
## Cd14 0.980 0.496 1.724797e-19
## Cebpb 0.987 0.556 1.161169e-16
## Pla2g7 0.709 0.103 1.042041e-14
## Pid1 0.788 0.231 3.501209e-14
## Ctsb 0.967 0.556 9.982281e-14
## C1qa 0.967 0.538 7.593105e-09table(macrophages.2.integrated$sample.id)##
## C0 C24
## 212 268colnames(conserved.markers.pos)## [1] "cluster" "gene" "max_pval" "minimump_p_val"
## [5] "marker.type" "C0_p_val" "C0_avg_log2FC" "C0_pct.1"
## [9] "C0_pct.2" "C0_p_val_adj" "C24_p_val" "C24_avg_log2FC"
## [13] "C24_pct.1" "C24_pct.2" "C24_p_val_adj"library(dplyr)
top1.pos <- conserved.markers.pos %>% group_by(cluster) %>% top_n(n = 1,
wt = -minimump_p_val)
top1.pos <- top1.pos %>% group_by(cluster) %>% top_n(n = 1, wt = -max_pval)
top1.pos <- top1.pos %>% group_by(cluster) %>% top_n(n = 1, wt = C24_avg_log2FC)
head(top1.pos)## # A tibble: 3 x 15
## # Groups: cluster [3]
## cluster gene max_pval minimump_p_val marker.type C0_p_val C0_avg_log2FC
## <chr> <chr> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 Macro1 Cd14 3.93e-13 2.41e-23 POS 3.93e-13 1.98
## 2 Macro2 Gpx3 7.94e- 8 3.02e-21 POS 7.94e- 8 1.24
## 3 Macro3 Pfkp 3.18e- 2 1.29e-31 POS 3.18e- 2 0.168
## # … with 8 more variables: C0_pct.1 <dbl>, C0_pct.2 <dbl>, C0_p_val_adj <dbl>,
## # C24_p_val <dbl>, C24_avg_log2FC <dbl>, C24_pct.1 <dbl>, C24_pct.2 <dbl>,
## # C24_p_val_adj <dbl>DefaultAssay(macrophages.2.integrated) <- "RNA"
f <- FeaturePlot(macrophages.2.integrated, features = top1.pos$gene,
min.cutoff = "q9", order = T, label = T, label.size = 8,
repel = T)
f
pdf(paste0(outdir, "/51_FeaturePlot_top_pos_markers_macrophages_2.pdf"),
width = 14, height = 14)
f
dev.off()## png
## 2DotPlot of top markers
DefaultAssay(macrophages.2.integrated) <- "RNA"
d <- DotPlot(object = macrophages.2.integrated, features = top1.pos$gene,
cols = c("#03539C", "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.2.subclusters") +
RotatedAxis()
d
pdf(paste0(outdir, "/52_DotPlot_top_cell_types_pos_markers_clusters.macrophages_2.pdf"),
width = 7, height = 7)
d
dev.off()## png
## 2DotPlot of top 10 markers
library(dplyr)
top10.pos <- conserved.markers.pos %>% group_by(cluster) %>%
top_n(n = 10, wt = -minimump_p_val)
top10.pos <- top10.pos %>% group_by(cluster) %>% top_n(n = 10,
wt = -max_pval)
top10.pos <- top10.pos %>% group_by(cluster) %>% top_n(n = 10,
wt = C24_avg_log2FC)
head(top10.pos)## # A tibble: 6 x 15
## # Groups: cluster [1]
## cluster gene max_pval minimump_p_val marker.type C0_p_val C0_avg_log2FC
## <chr> <chr> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 Macro1 Cd14 3.93e-13 2.41e-23 POS 3.93e-13 1.98
## 2 Macro1 Cebpb 1.07e- 2 1.62e-20 POS 1.07e- 2 0.668
## 3 Macro1 Pla2… 1.07e- 5 1.46e-18 POS 1.07e- 5 1.31
## 4 Macro1 Pid1 1.66e- 1 4.90e-18 POS 1.66e- 1 0.0460
## 5 Macro1 Ctsb 1.45e- 6 1.40e-17 POS 1.45e- 6 0.846
## 6 Macro1 C1qa 5.31e-13 1.48e-16 POS 7.40e-17 1.70
## # … with 8 more variables: C0_pct.1 <dbl>, C0_pct.2 <dbl>, C0_p_val_adj <dbl>,
## # C24_p_val <dbl>, C24_avg_log2FC <dbl>, C24_pct.1 <dbl>, C24_pct.2 <dbl>,
## # C24_p_val_adj <dbl>DefaultAssay(macrophages.2.integrated) <- "RNA"
d <- DotPlot(object = macrophages.2.integrated, features = top10.pos$gene,
cols = c("#03539C", "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.2.subclusters") +
RotatedAxis()
d
pdf(paste0(outdir, "/53_DotPlot_top10_cell_types_pos_markers_clusters.macrophages.2.pdf"),
width = 14, height = 7)
d
dev.off()## png
## 2Segregation of clusters by various sources of uninteresting variation
Explore additional metrics, such as the number of UMIs and genes per cell, and mitochondrial gene expression by UMAP.
# Determine metrics to plot present in
# tcells.integrated@meta.data
metrics <- c("nCount_RNA", "nFeature_RNA", "percent.mito")
f <- FeaturePlot(macrophages.2.integrated, reduction = "umap",
features = metrics, pt.size = 0.4, order = TRUE, min.cutoff = "q10",
label = TRUE)
f
pdf(paste0(outdir, "/54_FeaturePlot_macrophages.2.integrated_umap_metrics.pdf"),
width = 14, height = 14)
f
dev.off()## png
## 2Remove cells with mito expn > 10%
dim(macrophages.integrated.filtered)## [1] 14293 480keep_cell <- macrophages.integrated.filtered$percent.mito <=
10
summary(keep_cell)## Mode FALSE TRUE
## logical 132 348macrophages.integrated.filtered <- macrophages.integrated.filtered[,
keep_cell]
dim(macrophages.integrated.filtered)## [1] 14293 348New subclustering after filtering
DefaultAssay(macrophages.integrated.filtered) <- "RNA"
macrophages.3.list <- SplitObject(macrophages.integrated.filtered,
split.by = "sample.id")
macrophages.3.list## $C0
## An object of class Seurat
## 16293 features across 156 samples within 2 assays
## Active assay: RNA (14293 features, 0 variable features)
## 1 other assay present: integrated
## 2 dimensional reductions calculated: pca, umap
##
## $C24
## An object of class Seurat
## 16293 features across 192 samples within 2 assays
## Active assay: RNA (14293 features, 0 variable features)
## 1 other assay present: integrated
## 2 dimensional reductions calculated: pca, umapfor (i in 1:length(macrophages.3.list)) {
macrophages.3.list[[i]] <- FindVariableFeatures(macrophages.3.list[[i]],
selection.method = "vst", nfeatures = 2000, verbose = FALSE)
}k.filter <- min(200, min(sapply(macrophages.3.list, ncol)))
k.filter## [1] 156macrophages.3.anchors <- FindIntegrationAnchors(object.list = macrophages.3.list,
dims = 1:30, k.filter = k.filter)macrophages.3.integrated <- IntegrateData(anchorset = macrophages.3.anchors,
dims = 1:30)DefaultAssay(macrophages.3.integrated) <- "integrated"
macrophages.3.integrated <- ScaleData(macrophages.3.integrated,
verbose = T)dim(macrophages.3.integrated)## [1] 2000 348macrophages.3.integrated <- RunPCA(macrophages.3.integrated,
npcs = 50, verbose = T)system.time(macrophages.3.integrated <- JackStraw(macrophages.3.integrated,
num.replicate = 100, dims = 50))## user system elapsed
## 57.091 0.007 57.147macrophages.3.integrated <- ScoreJackStraw(macrophages.3.integrated,
dims = 1:50)
j <- JackStrawPlot(macrophages.3.integrated, dims = 1:50)
pdf(paste0(outdir, "/55_JackStrawPlot_macrophages_3.pdf"))
j
dev.off()## png
## 2macrophages.3.integrated <- FindNeighbors(macrophages.3.integrated,
dims = 1:42)
macrophages.3.integrated <- RunUMAP(macrophages.3.integrated,
dims = 1:42)for (i in seq(0, 2, 0.1)) {
macrophages.3.integrated <- FindClusters(macrophages.3.integrated,
resolution = i, verbose = F)
macrophages.3.integrated[[paste0("macrophages_3_subclusters_res.",
i)]] <- macrophages.3.integrated$seurat_clusters
}# install.packages('clustree')
library(clustree)
c <- clustree(macrophages.3.integrated, prefix = "macrophages_3_subclusters_res.")
c
pdf(paste0(outdir, "/56_clustree_macrophages_3.pdf"))
c
dev.off()## png
## 2Using the stability index to assess clusters The stability index from the SC3 package (Kiselev et al. 2017) measures the stability of clusters across resolutions and is automatically calculated when a clustering tree is built.
c <- clustree(macrophages.3.integrated, prefix = "macrophages_3_subclusters_res.",
node_colour = "sc3_stability")
pdf(paste0(outdir, "/57_clustree_macrophages_3_sc3_stability.pdf"))
c
dev.off()## png
## 2SOCS3, IL-1R2, IL1rn
Socs3, Il1r2, Il1rn
c <- clustree(macrophages.3.integrated, prefix = "macrophages_3_subclusters_res.",
node_colour = "Socs3", node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/58_clustree_macrophages_3_Socs3.pdf"))
c
dev.off()## png
## 2c <- clustree(macrophages.3.integrated, prefix = "macrophages_3_subclusters_res.",
node_colour = "Il1r2", node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/59_clustree_macrophages_3_Il1r2.pdf"))
c
dev.off()## png
## 2c <- clustree(macrophages.3.integrated, prefix = "macrophages_3_subclusters_res.",
node_colour = "Il1rn", node_colour_aggr = "mean")
c
pdf(paste0(outdir, "/60_clustree_macrophages_3_Il1rn.pdf"))
c
dev.off()## png
## 2RES 0.3
# install.packages('Polychrome')
library(Polychrome)
set.seed(5658807) # for reproducibility
levels <- levels(as.factor(macrophages.3.integrated$macrophages_3_subclusters_res.0.3))
rainbow2.6c <- c("#03539C", "#12999E", "#B7CE05", "#FAEB09",
"#FDA908", "#E82564")
colors.macrophages.3.subclusters <- createPalette(length(levels),
rainbow2.6c, M = 1000)
names(colors.macrophages.3.subclusters) <- levels
colors.macrophages.3.subclusters## 0 1 2 <NA> <NA> <NA> <NA> <NA>
## "#1662B5" "#0D9A9D" "#C2DB1C" "#FBE22E" "#FCA700" "#F33877" "#FA35F5" "#FFB1DB"
## <NA> <NA> <NA>
## "#40F897" "#63691C" "#4F00FF"slices <- rep(1, length(colors.macrophages.3.subclusters))
pie(slices, col = colors.macrophages.3.subclusters, labels = names(colors.macrophages.3.subclusters))
Idents(macrophages.3.integrated) <- "macrophages_3_subclusters_res.0.3"
DimPlot(macrophages.3.integrated, reduction = "umap", pt.size = 1,
label = T, label.size = 8, cols = colors.macrophages.3.subclusters)
pdf(paste0(outdir, "/61_DimPlot_umap_macrophages_3_subclusters_res.0.3.pdf"))
DimPlot(macrophages.3.integrated, reduction = "umap", pt.size = 1,
label = T, label.size = 6, cols = colors.macrophages.3.subclusters)
dev.off()## png
## 2colors.samples <- c("#12999E", "#FDA908")
names(colors.samples) <- c("C0", "C24")
macrophages.3.integrated$sample.id <- factor(macrophages.3.integrated$sample.id,
levels = names(colors.samples))
p <- DimPlot(macrophages.3.integrated, reduction = "umap", group.by = "sample.id",
cols = colors.samples, label = F)
pdf(paste0(outdir, "/62_UMAP_macrophages_3_samples.pdf"), width = 10,
height = 9)
p
dev.off()## png
## 2pdf(paste0(outdir, "/63_UMAP_macrophages_3_samples_noLegend.pdf"),
width = 7, height = 7)
p + NoLegend()
dev.off()## png
## 2Rename subclusters
new.cluster.ids <- paste0("Macro", 1:length(levels(as.factor(macrophages.3.integrated$macrophages_3_subclusters_res.0.3))))
names(new.cluster.ids) <- levels(macrophages.3.integrated)
macrophages.3.integrated <- RenameIdents(macrophages.3.integrated,
new.cluster.ids)
macrophages.3.integrated[["renamed.macrophages.3.subclusters"]] <- macrophages.3.integrated@active.identIdents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters"
names(colors.macrophages.3.subclusters) <- levels(macrophages.3.integrated$renamed.macrophages.3.subclusters)
DimPlot(macrophages.3.integrated, reduction = "umap", pt.size = 1,
label = T, label.size = 8, cols = colors.macrophages.3.subclusters)
pdf(paste0(outdir, "/64_DimPlot_umap_renamed_macrophages_3_subclusters.pdf"))
DimPlot(macrophages.3.integrated, reduction = "umap", pt.size = 1,
label = T, label.size = 6, cols = colors.macrophages.3.subclusters)
dev.off()## png
## 2Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters"
DimPlot(macrophages.3.integrated, reduction = "umap", pt.size = 1,
label = F, label.size = 8, cols = colors.macrophages.3.subclusters) +
NoLegend()
pdf(paste0(outdir, "/65_DimPlot_umap_renamed_macrophages_3_subclusters_noLabels_noLegend.pdf"))
DimPlot(macrophages.3.integrated, reduction = "umap", pt.size = 1,
label = F, label.size = 6, cols = colors.macrophages.3.subclusters) +
NoLegend()
dev.off()## png
## 2saveRDS(macrophages.3.integrated, paste0(outdir, "/66.macrophages.3.integrated.rds"))Conserved markers
t <- table(macrophages.3.integrated$renamed.macrophages.3.subclusters,
macrophages.3.integrated$sample.id)
t##
## C0 C24
## Macro1 108 118
## Macro2 34 42
## Macro3 14 32openxlsx::write.xlsx(as.data.frame.matrix(t), file = paste0(outdir,
"/67_Sample_macrophages_3_subcluster_composition.xlsx"),
rowNames = T, colNames = T)DefaultAssay(macrophages.3.integrated) <- "RNA"
Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters"
subclusters <- levels(as.factor(macrophages.3.integrated$renamed.macrophages.3.subclusters))
conserved.markers <- data.frame(matrix(ncol = 14))
for (c in subclusters) {
print(c)
markers.c <- FindConservedMarkers(macrophages.3.integrated,
ident.1 = c, grouping.var = "sample.id", verbose = T,
logfc.threshold = -Inf, min.pct = -Inf, min.cells.feature = -Inf,
min.cells.group = -Inf)
markers.c <- cbind(data.frame(cluster = rep(c, dim(markers.c)[1]),
gene = rownames(markers.c)), markers.c)
write.table(markers.c, file = paste0(outdir, "/68_markers_",
c, "_macrophages_3_subclusters.txt"))
colnames(conserved.markers) <- colnames(markers.c)
conserved.markers <- rbind(conserved.markers, markers.c)
}## [1] "Macro1"
## [1] "Macro2"
## [1] "Macro3"# for (c in 1:4) { subcluster <- paste0('AM', c) markers.c <-
# read.table(file = paste0(outdir, '/31_markers_',
# subcluster, '_AM_subclusters.txt'))
# colnames(conserved.markers) <- colnames(markers.c)
# conserved.markers <- rbind(conserved.markers, markers.c) }
# head(conserved.markers)# markers.AM5 <- read.table(file = paste0(outdir,
# '/31_markers_', 'AM5', '_AM_subclusters.txt'))
# colnames(markers.AM5) colnames(conserved.markers)# for ( i in 1:10 ) { markers.AM5[, 34+i] <- rep(NA,
# dim(markers.AM5)[1]) } markers.AM5 <- markers.AM5[, c(1:12,
# 35:44, 13:34)] colnames(markers.AM5) <-
# colnames(conserved.markers) head(markers.AM5)# conserved.markers <- rbind(conserved.markers, markers.AM5)conserved.markers <- conserved.markers[-1, ]
conserved.markers <- conserved.markers[, c(1:2, 13:14, 3:12)]
head(conserved.markers)## cluster gene max_pval minimump_p_val C0_p_val C0_avg_log2FC
## Cd14 Macro1 Cd14 3.903045e-10 2.835224e-21 3.903045e-10 1.9896992
## C1qa Macro1 C1qa 7.937402e-15 9.276577e-17 4.638289e-17 2.2670625
## C1qc Macro1 C1qc 7.233602e-13 3.884857e-15 7.233602e-13 1.6955969
## Cebpb Macro1 Cebpb 3.638433e-01 7.452346e-15 3.638433e-01 0.1487639
## Sdc4 Macro1 Sdc4 4.119571e-02 1.039183e-14 4.119571e-02 0.3247144
## C1qb Macro1 C1qb 8.047736e-12 2.132823e-14 1.066411e-14 1.9546166
## C0_pct.1 C0_pct.2 C0_p_val_adj C24_p_val C24_avg_log2FC C24_pct.1
## Cd14 0.796 0.271 5.578622e-06 1.417612e-21 2.201101 1.000
## C1qa 0.972 0.396 6.629506e-13 7.937402e-15 1.509209 0.992
## C1qc 0.981 0.375 1.033899e-08 1.942429e-15 1.617483 0.975
## Cebpb 0.185 0.125 1.000000e+00 3.726173e-15 1.545982 0.983
## Sdc4 0.306 0.146 1.000000e+00 5.195913e-15 1.213547 1.000
## C1qb 0.972 0.417 1.524222e-10 8.047736e-12 1.113736 1.000
## C24_pct.2 C24_p_val_adj
## Cd14 0.541 2.026193e-17
## C1qa 0.405 1.134493e-10
## C1qc 0.365 2.776313e-11
## Cebpb 0.581 5.325819e-11
## Sdc4 0.784 7.426519e-11
## C1qb 0.486 1.150263e-07Only top markers that are positive markers
conserved.markers$marker.type <- apply(conserved.markers, 1, function(x) {
y <- as.numeric(x)
if ( (if (is.na(y[6])) {TRUE} else {y[6]>0})
& (if (is.na(y[11])) {TRUE} else {y[11]>0})
# & (if (is.na(y[16])) {TRUE} else {y[16]>0})
# & (if (is.na(y[21])) {TRUE} else {y[21]>0})
# & (if (is.na(y[26])) {TRUE} else {y[26]>0})
# & (if (is.na(y[31])) {TRUE} else {y[31]>0})
# & (if (is.na(y[36])) {TRUE} else {y[36]>0})
# & (if (is.na(y[41])) {TRUE} else {y[41]>0})
)
{"POS"}
else if ( ( if (is.na(y[6])) {TRUE} else {y[6]<0})
& (if (is.na(y[11])) {TRUE} else {y[11]<0})
# & (if (is.na(y[16])) {TRUE} else {y[16]<0})
# & (if (is.na(y[21])) {TRUE} else {y[21]<0})
# & (if (is.na(y[26])) {TRUE} else {y[26]<0})
# & (if (is.na(y[31])) {TRUE} else {y[31]<0})
# & (if (is.na(y[36])) {TRUE} else {y[36]<0})
# & (if (is.na(y[41])) {TRUE} else {y[41]<0})
)
{"NEG"}
else {"UND"}
})
conserved.markers <- conserved.markers[, c(1:4, 15, 5:14)]openxlsx::write.xlsx(conserved.markers, paste0(outdir, "/69_conserved_markers_macrophages_3.xlsx"),
colNames = T)
head(conserved.markers)## cluster gene max_pval minimump_p_val marker.type C0_p_val
## Cd14 Macro1 Cd14 3.903045e-10 2.835224e-21 POS 3.903045e-10
## C1qa Macro1 C1qa 7.937402e-15 9.276577e-17 POS 4.638289e-17
## C1qc Macro1 C1qc 7.233602e-13 3.884857e-15 POS 7.233602e-13
## Cebpb Macro1 Cebpb 3.638433e-01 7.452346e-15 POS 3.638433e-01
## Sdc4 Macro1 Sdc4 4.119571e-02 1.039183e-14 POS 4.119571e-02
## C1qb Macro1 C1qb 8.047736e-12 2.132823e-14 POS 1.066411e-14
## C0_avg_log2FC C0_pct.1 C0_pct.2 C0_p_val_adj C24_p_val C24_avg_log2FC
## Cd14 1.9896992 0.796 0.271 5.578622e-06 1.417612e-21 2.201101
## C1qa 2.2670625 0.972 0.396 6.629506e-13 7.937402e-15 1.509209
## C1qc 1.6955969 0.981 0.375 1.033899e-08 1.942429e-15 1.617483
## Cebpb 0.1487639 0.185 0.125 1.000000e+00 3.726173e-15 1.545982
## Sdc4 0.3247144 0.306 0.146 1.000000e+00 5.195913e-15 1.213547
## C1qb 1.9546166 0.972 0.417 1.524222e-10 8.047736e-12 1.113736
## C24_pct.1 C24_pct.2 C24_p_val_adj
## Cd14 1.000 0.541 2.026193e-17
## C1qa 0.992 0.405 1.134493e-10
## C1qc 0.975 0.365 2.776313e-11
## Cebpb 0.983 0.581 5.325819e-11
## Sdc4 1.000 0.784 7.426519e-11
## C1qb 1.000 0.486 1.150263e-07saveRDS(macrophages.3.integrated, paste0(outdir, "/70.macrophages.3.integrated.rds"))conserved.markers.pos <- conserved.markers[conserved.markers$marker.type ==
"POS", ]
head(conserved.markers.pos)## cluster gene max_pval minimump_p_val marker.type C0_p_val
## Cd14 Macro1 Cd14 3.903045e-10 2.835224e-21 POS 3.903045e-10
## C1qa Macro1 C1qa 7.937402e-15 9.276577e-17 POS 4.638289e-17
## C1qc Macro1 C1qc 7.233602e-13 3.884857e-15 POS 7.233602e-13
## Cebpb Macro1 Cebpb 3.638433e-01 7.452346e-15 POS 3.638433e-01
## Sdc4 Macro1 Sdc4 4.119571e-02 1.039183e-14 POS 4.119571e-02
## C1qb Macro1 C1qb 8.047736e-12 2.132823e-14 POS 1.066411e-14
## C0_avg_log2FC C0_pct.1 C0_pct.2 C0_p_val_adj C24_p_val C24_avg_log2FC
## Cd14 1.9896992 0.796 0.271 5.578622e-06 1.417612e-21 2.201101
## C1qa 2.2670625 0.972 0.396 6.629506e-13 7.937402e-15 1.509209
## C1qc 1.6955969 0.981 0.375 1.033899e-08 1.942429e-15 1.617483
## Cebpb 0.1487639 0.185 0.125 1.000000e+00 3.726173e-15 1.545982
## Sdc4 0.3247144 0.306 0.146 1.000000e+00 5.195913e-15 1.213547
## C1qb 1.9546166 0.972 0.417 1.524222e-10 8.047736e-12 1.113736
## C24_pct.1 C24_pct.2 C24_p_val_adj
## Cd14 1.000 0.541 2.026193e-17
## C1qa 0.992 0.405 1.134493e-10
## C1qc 0.975 0.365 2.776313e-11
## Cebpb 0.983 0.581 5.325819e-11
## Sdc4 1.000 0.784 7.426519e-11
## C1qb 1.000 0.486 1.150263e-07table(macrophages.3.integrated$sample.id)##
## C0 C24
## 156 192colnames(conserved.markers.pos)## [1] "cluster" "gene" "max_pval" "minimump_p_val"
## [5] "marker.type" "C0_p_val" "C0_avg_log2FC" "C0_pct.1"
## [9] "C0_pct.2" "C0_p_val_adj" "C24_p_val" "C24_avg_log2FC"
## [13] "C24_pct.1" "C24_pct.2" "C24_p_val_adj"library(dplyr)
top1.pos <- conserved.markers.pos %>% group_by(cluster) %>% top_n(n = 1,
wt = -minimump_p_val)
top1.pos <- top1.pos %>% group_by(cluster) %>% top_n(n = 1, wt = -max_pval)
top1.pos <- top1.pos %>% group_by(cluster) %>% top_n(n = 1, wt = C24_avg_log2FC)
head(top1.pos)## # A tibble: 3 x 15
## # Groups: cluster [3]
## cluster gene max_pval minimump_p_val marker.type C0_p_val C0_avg_log2FC
## <chr> <chr> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 Macro1 Cd14 3.90e-10 2.84e-21 POS 3.90e-10 1.99
## 2 Macro2 Pfkp 1.52e- 2 2.28e-23 POS 1.52e- 2 0.0989
## 3 Macro3 Kap 4.60e- 6 7.21e-11 POS 4.60e- 6 1.77
## # … with 8 more variables: C0_pct.1 <dbl>, C0_pct.2 <dbl>, C0_p_val_adj <dbl>,
## # C24_p_val <dbl>, C24_avg_log2FC <dbl>, C24_pct.1 <dbl>, C24_pct.2 <dbl>,
## # C24_p_val_adj <dbl>DefaultAssay(macrophages.3.integrated) <- "RNA"
f <- FeaturePlot(macrophages.3.integrated, features = top1.pos$gene,
min.cutoff = "q9", order = T, label = T, label.size = 8,
repel = T)
f
pdf(paste0(outdir, "/71_FeaturePlot_top_pos_markers_macrophages_3.pdf"),
width = 14, height = 14)
f
dev.off()## png
## 2DotPlot of top markers
DefaultAssay(macrophages.3.integrated) <- "RNA"
d <- DotPlot(object = macrophages.3.integrated, features = top1.pos$gene,
cols = c("#03539C", "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.3.subclusters") +
RotatedAxis()
d
pdf(paste0(outdir, "/72_DotPlot_top_cell_types_pos_markers_clusters.macrophages_3.pdf"),
width = 7, height = 7)
d
dev.off()## png
## 2DotPlot of top 10 markers
library(dplyr)
top10.pos <- conserved.markers.pos %>% group_by(cluster) %>%
top_n(n = 10, wt = -minimump_p_val)
top10.pos <- top10.pos %>% group_by(cluster) %>% top_n(n = 10,
wt = -max_pval)
top10.pos <- top10.pos %>% group_by(cluster) %>% top_n(n = 10,
wt = C24_avg_log2FC)
head(top10.pos)## # A tibble: 6 x 15
## # Groups: cluster [1]
## cluster gene max_pval minimump_p_val marker.type C0_p_val C0_avg_log2FC
## <chr> <chr> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 Macro1 Cd14 3.90e-10 2.84e-21 POS 3.90e-10 1.99
## 2 Macro1 C1qa 7.94e-15 9.28e-17 POS 4.64e-17 2.27
## 3 Macro1 C1qc 7.23e-13 3.88e-15 POS 7.23e-13 1.70
## 4 Macro1 Cebpb 3.64e- 1 7.45e-15 POS 3.64e- 1 0.149
## 5 Macro1 Sdc4 4.12e- 2 1.04e-14 POS 4.12e- 2 0.325
## 6 Macro1 C1qb 8.05e-12 2.13e-14 POS 1.07e-14 1.95
## # … with 8 more variables: C0_pct.1 <dbl>, C0_pct.2 <dbl>, C0_p_val_adj <dbl>,
## # C24_p_val <dbl>, C24_avg_log2FC <dbl>, C24_pct.1 <dbl>, C24_pct.2 <dbl>,
## # C24_p_val_adj <dbl>DefaultAssay(macrophages.3.integrated) <- "RNA"
d <- DotPlot(object = macrophages.3.integrated, features = top10.pos$gene,
cols = c("#03539C", "#E82564"), dot.scale = 8, group.by = "renamed.macrophages.3.subclusters") +
RotatedAxis()
d
pdf(paste0(outdir, "/73_DotPlot_top10_cell_types_pos_markers_clusters.macrophages.3.pdf"),
width = 14, height = 7)
d
dev.off()## png
## 2Segregation of clusters by various sources of uninteresting variation
Explore additional metrics, such as the number of UMIs and genes per cell, and mitochondrial gene expression by UMAP.
# Determine metrics to plot present in
# tcells.integrated@meta.data
metrics <- c("nCount_RNA", "nFeature_RNA", "percent.mito")
f <- FeaturePlot(macrophages.3.integrated, reduction = "umap",
features = metrics, pt.size = 0.4, order = TRUE, min.cutoff = "q10",
label = TRUE)
f
pdf(paste0(outdir, "/74_FeaturePlot_macrophages.3.integrated_umap_metrics.pdf"),
width = 14, height = 14)
f
dev.off()## png
## 2Heatmap
macrophages.3.integrated <- ScaleData(macrophages.3.integrated,
verbose = T, features = rownames(macrophages.3.integrated))
d <- DoHeatmap(macrophages.3.integrated, features = top10.pos$gene,
group.by = "renamed.macrophages.3.subclusters", group.colors = colors.macrophages.3.subclusters) +
NoLegend()
d
pdf(paste0(outdir, "/75_DoHeatmap_top10_pos_macrophages_3.pdf"),
width = 15, height = 13)
d
dev.off()## png
## 2DotPlot with known/useful markers
immune.cells <- c("Ptprc", "Spi1", "Cd3g")
b.cells <- c("Cd19", "Cd79a", "Cd79b")
t.cells <- "Cd3e"
nk.cells <- c("Ncr1")
myeloid.cells <- "Cd68"
neutrophils <- c("Csf3r", "S100a8", "S100a9", "Retnlg")
monocytes.macrophages.dendritic.cells <- "Csf1r"
dc <- c("Siglech", "Cd209a")
dc1 <- "Xcr1"
dc2 <- c("Ccl18", "Nr4a3", "Tnfsf12")
dc3 <- "Fscn1"
pdc <- "Tcf4"
monocytes <- c("Ly6c1", "Ly6c2", "Ccr2", "Cx3cr1", "Fcgr4")
macrophages <- c("Mrc1", "Lyz2", "Adgre1", "Axl", "Cxcl2", "Trem2",
"C1qc", "Fcgr1", "C5ar1")
am.markers <- c("Pparg", "Itgax", "Ear1", "Ear2", "Ear3", "Car4",
"Siglec5", "SiglecF")
im.markers <- c("Itgam", "Cx3cr1", "Ccl12", "Bhlme40", "Trem2",
"Cxcl12", "Csf1")
peribronchial.macrophages <- c("Ccr2", "Bhlhe40", "Bhlhe41")
perivascular.macrophages <- c("Lyve1", "F13a1")
proliferating.cells <- "Mki67"
mesenchymal.cells <- c("Col1a1", "Col1a2", "Pdgfra")
epithelial.cells <- c("Epcam", "Cdh1")
mast.cells <- "Mcpt8"
endothelial.cells <- c("Pecam1", "Emcn")markers <- unique(c(immune.cells, b.cells, t.cells, nk.cells,
myeloid.cells, neutrophils, monocytes.macrophages.dendritic.cells,
dc, dc1, dc2, dc3, pdc, monocytes, macrophages, am.markers,
im.markers, peribronchial.macrophages, perivascular.macrophages,
proliferating.cells, mesenchymal.cells, epithelial.cells,
mast.cells, endothelial.cells))Dot plot
DefaultAssay(macrophages.3.integrated) <- "RNA"
d <- DotPlot(object = macrophages.3.integrated, features = markers,
cols = c("#03539C", "#E82564"), dot.scale = 15, group.by = "renamed.macrophages.3.subclusters") +
RotatedAxis()
d
pdf(paste0(outdir, "/76_DotPlot_known_markers_macrophages_3.pdf"),
width = 18, height = 4)
d
dev.off()## png
## 2Expression of genes of interest: Socs3, Il1r2, Il1rn
goi <- c("Socs3", "Il1r2", "Il1rn")DefaultAssay(macrophages.3.integrated) <- "RNA"
f <- FeaturePlot(macrophages.3.integrated, features = goi, min.cutoff = "q9",
order = T, label = T, label.size = 8, repel = T)
f
pdf(paste0(outdir, "/77_FeaturePlot_goi_macrophages_3.pdf"),
width = 14, height = 14)
f
dev.off()DefaultAssay(macrophages.3.integrated) <- "RNA"
d <- DotPlot(object = macrophages.3.integrated, features = goi,
cols = c("#03539C", "#E82564"), dot.scale = 15, group.by = "renamed.macrophages.3.subclusters") +
RotatedAxis()
d
pdf(paste0(outdir, "/78_DotPlot_goi_macrophages.pdf"), width = 6,
height = 4)
d
dev.off()## png
## 2DIFFERENTIAL EXPRESSION ANALYSES
differential expression studies for the following groups: 1. C24 vs C0 within each cell type (from manual annotation annotation.1)
1. C24 vs C0 within each cell type (from manual annotation annotation.1)
macrophages.3.integrated$renamed.macrophages.3.subclusters.sample.id <- paste(macrophages.3.integrated$renamed.macrophages.3.subclusters,
macrophages.3.integrated$sample.id, sep = "_")
levels(factor(macrophages.3.integrated$renamed.macrophages.3.subclusters.sample.id))## [1] "Macro1_C0" "Macro1_C24" "Macro2_C0" "Macro2_C24" "Macro3_C0"
## [6] "Macro3_C24"DefaultAssay(macrophages.3.integrated) <- "RNA"
output <- paste0(outdir, "/79_DE_macrophages_3_subclusters_C24_vs_C0_")
levels(factor(macrophages.3.integrated$renamed.macrophages.3.subclusters.sample.id))## [1] "Macro1_C0" "Macro1_C24" "Macro2_C0" "Macro2_C24" "Macro3_C0"
## [6] "Macro3_C24"for (cell.type in levels(factor(macrophages.3.integrated$renamed.macrophages.3.subclusters))) {
try({
print(cell.type)
ident1 <- paste0(cell.type, "_C24")
ident2 <- paste0(cell.type, "_C0")
Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters.sample.id"
condition.diffgenes <- FindMarkers(macrophages.3.integrated,
ident.1 = ident1, ident.2 = ident2, logfc.threshold = -Inf,
test.use = "wilcox", min.pct = -Inf, verbose = T,
min.cells.feature = -Inf, min.cells.group = -Inf)
condition.diffgenes.adjpval.0.05 <- condition.diffgenes[condition.diffgenes$p_val_adj <
0.05, ]
condition.diffgenes.adjpval.0.05.upregulated.sorted <- condition.diffgenes.adjpval.0.05[condition.diffgenes.adjpval.0.05$avg_log2FC >
0, ]
condition.diffgenes.adjpval.0.05.upregulated.sorted <- condition.diffgenes.adjpval.0.05.upregulated.sorted[order(condition.diffgenes.adjpval.0.05.upregulated.sorted$avg_log2FC,
decreasing = T), ]
condition.diffgenes.adjpval.0.05.downregulated.sorted <- condition.diffgenes.adjpval.0.05[condition.diffgenes.adjpval.0.05$avg_log2FC <
0, ]
condition.diffgenes.adjpval.0.05.downregulated.sorted <- condition.diffgenes.adjpval.0.05.downregulated.sorted[order(condition.diffgenes.adjpval.0.05.downregulated.sorted$avg_log2FC,
decreasing = F), ]
list_of_datasets <- list(sig.upregulated.genes.ordered.by.avg.logFC = condition.diffgenes.adjpval.0.05.upregulated.sorted,
sig.downregulated.genes.ordered.by.avg.logFC = condition.diffgenes.adjpval.0.05.downregulated.sorted,
all.genes.ordered.by.adj.pval = condition.diffgenes)
openxlsx::write.xlsx(list_of_datasets, file = paste0(output,
cell.type, ".xlsx"), row.names = T)
})
}## [1] "Macro1"
## [1] "Macro2"
## [1] "Macro3"for (cell.type in levels(factor(macrophages.3.integrated$renamed.macrophages.3.subclusters))) {
try({
print(cell.type)
filename <- paste0(output, cell.type, ".xlsx")
sheets <- openxlsx::getSheetNames(filename)
SheetList <- lapply(sheets, openxlsx::read.xlsx, xlsxFile = filename)
names(SheetList) <- sheets
condition.diffgenes.adjpval.0.05.upregulated.sorted <- SheetList[[1]]
rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted) <- condition.diffgenes.adjpval.0.05.upregulated.sorted[,
1]
condition.diffgenes.adjpval.0.05.upregulated.sorted <- condition.diffgenes.adjpval.0.05.upregulated.sorted[,
-1]
condition.diffgenes.adjpval.0.05.downregulated.sorted <- SheetList[[2]]
rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted) <- condition.diffgenes.adjpval.0.05.downregulated.sorted[,
1]
condition.diffgenes.adjpval.0.05.downregulated.sorted <- condition.diffgenes.adjpval.0.05.downregulated.sorted[,
-1]
Idents(annotated) <- "renamed.macrophages.3.subclusters"
try({
f <- FeaturePlot(macrophages.3.integrated, features = rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted)[1],
split.by = "sample.id", cols = c("grey", "#E82564"),
order = T, min.cutoff = "q10", max.cutoff = "q90")
pdf(paste0(output, cell.type, "_FeaturePlot_top_upregulated_split_by_sample.pdf"),
width = 14, height = 7)
print(f)
dev.off()
})
try({
f <- FeaturePlot(macrophages.3.integrated, features = rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted)[1],
split.by = "sample.id", cols = c("grey", "#E82564"),
order = T, min.cutoff = "q10", max.cutoff = "q90")
pdf(paste0(output, cell.type, "_FeaturePlot_top_downregulated_split_by_sample.pdf"),
width = 14, height = 7)
print(f)
dev.off()
})
try({
plot <- VlnPlot(macrophages.3.integrated, features = rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted)[1],
group.by = "renamed.macrophages.3.subclusters",
split.by = "sample.id", pt.size = 0, combine = FALSE,
cols = colors.samples)
pdf(paste0(output, cell.type, "_VlnPlot_top_upregulated_group_by_cluster_split_by_sample.pdf"),
width = 11, height = 8.5)
print(plot)
dev.off()
})
try({
plot <- VlnPlot(macrophages.3.integrated, features = rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted)[1],
group.by = "renamed.macrophages.3.subclusters",
split.by = "sample.id", pt.size = 0, combine = FALSE,
cols = colors.samples)
pdf(paste0(output, cell.type, "_VlnPlot_top_downregulated_group_by_cluster_split_by_sample.pdf"),
width = 11, height = 8.5)
print(plot)
dev.off()
})
## DotPlot with 50 most upregulated genes
Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters.sample.id"
try({
d <- DotPlot(macrophages.3.integrated, features = rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted)[1:50],
cols = c("#03539C", "#E82564"), dot.scale = 8,
group.by = "renamed.macrophages.3.subclusters.sample.id") +
RotatedAxis()
pdf(paste0(output, cell.type, "_DotPlot_top50_upregulated_groupby_annotation.1.sample.id.pdf"),
width = (4 + min(50, length(rownames(condition.diffgenes.adjpval.0.05.upregulated.sorted))) *
15/50), height = 17)
print(d)
dev.off()
})
## DotPlot with 50 most downregulated genes
Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters.sample.id"
try({
d <- DotPlot(macrophages.3.integrated, features = rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted)[1:50],
cols = c("#03539C", "#E82564"), dot.scale = 8,
group.by = "renamed.macrophages.3.subclusters.sample.id") +
RotatedAxis()
pdf(paste0(output, cell.type, "_DotPlot_top50_downregulated_groupby_annotation.1.sample.id.pdf"),
width = (4 + min(50, length(rownames(condition.diffgenes.adjpval.0.05.downregulated.sorted))) *
15/50), height = 17)
print(d)
dev.off()
})
})
}## [1] "Macro1"## [[1]]##
## [[1]]## [1] "Macro2"## [[1]]##
## [[1]]## [1] "Macro3"## Error : None of the requested features were found: NA in slot data
## [[1]]##
## Error in FetchData(object = object, vars = features, slot = slot) :
## None of the requested variables were found: NA## Error in h(simpleError(msg, call)) :
## error in evaluating the argument 'x' in selecting a method for function 'mean': non-numeric argument to mathematical functionPlots for genes of interest
for (gene in goi) {
try({
Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters"
try({
f <- FeaturePlot(macrophages.3.integrated, features = gene,
split.by = "sample.id", cols = c("grey", "#E82564"),
order = T, min.cutoff = "q10", max.cutoff = "q90",
label = T)
pdf(paste0(outdir, "/80_", gene, "_FeaturePlot_split_by_sample.pdf"),
width = 14, height = 7)
print(f)
dev.off()
})
try({
plot <- VlnPlot(macrophages.3.integrated, features = gene,
group.by = "renamed.macrophages.3.subclusters",
split.by = "sample.id", pt.size = 0, combine = FALSE,
cols = colors.samples)
pdf(paste0(outdir, "/81_", gene, "_VlnPlot_group_by_cluster_split_by_sample.pdf"),
width = 11, height = 8.5)
print(plot)
dev.off()
})
})
}## [[1]]## [[1]]## [[1]]Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters.sample.id"
try({
d <- DotPlot(macrophages.3.integrated, features = goi, cols = c("#03539C",
"#E82564"), dot.scale = 8, group.by = "renamed.macrophages.3.subclusters.sample.id") +
RotatedAxis()
pdf(paste0(outdir, "/82_", paste(goi, sep = "_"), "_DotPlot_groupby_annotation.1.sample.id.pdf"),
width = (4 + length(goi) * 15/50), height = 17)
print(d)
dev.off()
})## png
## 2Plots for more genes of interest: F480 (antibody), CD11b, CD11c, CD64, MHCII
In mouse: F4/80 (antibody) > Emr1, Adgre1 CD11b -> Itgam CD11c -> Itgax CD64 -> Fcgr1 MHCII -> H2-Ab
goi <- c("Adgre1", "Itgam", "Itgax", "Fcgr1", "H2-Ab1")DefaultAssay(macrophages.3.integrated) <- "RNA"
for (gene in goi) {
c <- clustree(macrophages.3.integrated, prefix = "macrophages_3_subclusters_res.",
node_colour = gene, node_colour_aggr = "mean")
pdf(paste0(outdir, "/83_clustree_", gene, "_macrophages.pdf"))
print(c)
dev.off()
}for (gene in goi) {
try({
Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters"
try({
f <- FeaturePlot(macrophages.3.integrated, features = gene,
split.by = "sample.id", cols = c("grey", "#E82564"),
order = T, min.cutoff = "q10", max.cutoff = "q90",
label = T)
pdf(paste0(outdir, "/84_", gene, "_FeaturePlot_split_by_sample.pdf"),
width = 14, height = 7)
print(f)
dev.off()
})
try({
plot <- VlnPlot(macrophages.3.integrated, features = gene,
group.by = "renamed.macrophages.3.subclusters",
split.by = "sample.id", pt.size = 0, combine = FALSE,
cols = colors.samples)
pdf(paste0(outdir, "/85_", gene, "_VlnPlot_group_by_cluster_split_by_sample.pdf"),
width = 11, height = 8.5)
print(plot)
dev.off()
})
})
}## [[1]]## [[1]]## [[1]]## [[1]]## [[1]]Idents(macrophages.3.integrated) <- "renamed.macrophages.3.subclusters.sample.id"
try({
d <- DotPlot(macrophages.3.integrated, features = goi, cols = c("#03539C",
"#E82564"), dot.scale = 8, group.by = "renamed.macrophages.3.subclusters.sample.id") +
RotatedAxis()
pdf(paste0(outdir, "/86_", paste(goi, sep = "_"), "_DotPlot_groupby_macrophages_subclusters.sample.id.pdf"),
width = (4 + length(goi) * 15/50), height = 17)
print(d)
dev.off()
})## png
## 2Session Information
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.8 (Santiago)
##
## Matrix products: default
## BLAS: /gpfs/fs1/data/omicscore/Privratsky-Privratsky-20210215/scripts/conda/envs/privratsky/lib/libblas.so.3.8.0
## LAPACK: /gpfs/fs1/data/omicscore/Privratsky-Privratsky-20210215/scripts/conda/envs/privratsky/lib/liblapack.so.3.8.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] Polychrome_1.2.6 clustree_0.4.3
## [3] ggraph_2.0.4 scater_1.18.0
## [5] SingleCellExperiment_1.12.0 data.table_1.13.6
## [7] gridExtra_2.3 forcats_0.5.1
## [9] stringr_1.4.0 purrr_0.3.4
## [11] readr_1.4.0 tidyr_1.1.2
## [13] tibble_3.0.6 tidyverse_1.3.0
## [15] SingleR_1.4.0 celldex_1.0.0
## [17] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [19] GenomicRanges_1.42.0 GenomeInfoDb_1.26.0
## [21] IRanges_2.24.0 S4Vectors_0.28.0
## [23] BiocGenerics_0.36.0 MatrixGenerics_1.2.0
## [25] matrixStats_0.58.0 ggsci_2.9
## [27] cowplot_1.1.1 plyr_1.8.6
## [29] patchwork_1.1.1 dplyr_1.0.4
## [31] ggplot2_3.3.3 SeuratObject_4.0.0
## [33] Seurat_4.0.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.1.4 reticulate_1.18
## [3] tidyselect_1.1.0 RSQLite_2.2.4
## [5] AnnotationDbi_1.52.0 htmlwidgets_1.5.3
## [7] grid_4.0.3 BiocParallel_1.24.0
## [9] Rtsne_0.15 munsell_0.5.0
## [11] codetools_0.2-18 ica_1.0-2
## [13] future_1.21.0 miniUI_0.1.1.1
## [15] withr_2.4.1 colorspace_2.0-0
## [17] highr_0.8 knitr_1.31
## [19] rstudioapi_0.13 ROCR_1.0-11
## [21] tensor_1.5 gbRd_0.4-11
## [23] listenv_0.8.0 Rdpack_2.1
## [25] labeling_0.4.2 GenomeInfoDbData_1.2.4
## [27] polyclip_1.10-0 farver_2.0.3
## [29] bit64_4.0.5 parallelly_1.23.0
## [31] vctrs_0.3.6 generics_0.1.0
## [33] xfun_0.20 BiocFileCache_1.14.0
## [35] R6_2.5.0 graphlayouts_0.7.1
## [37] ggbeeswarm_0.6.0 rsvd_1.0.3
## [39] bitops_1.0-6 spatstat.utils_2.0-0
## [41] cachem_1.0.4 DelayedArray_0.16.0
## [43] assertthat_0.2.1 promises_1.2.0.1
## [45] scales_1.1.1 beeswarm_0.2.3
## [47] gtable_0.3.0 beachmat_2.6.0
## [49] globals_0.14.0 goftest_1.2-2
## [51] tidygraph_1.2.0 rlang_0.4.10
## [53] scatterplot3d_0.3-41 splines_4.0.3
## [55] lazyeval_0.2.2 checkmate_2.0.0
## [57] broom_0.7.5 BiocManager_1.30.10
## [59] yaml_2.2.1 reshape2_1.4.4
## [61] abind_1.4-5 modelr_0.1.8
## [63] backports_1.2.1 httpuv_1.5.5
## [65] tools_4.0.3 ellipsis_0.3.1
## [67] RColorBrewer_1.1-2 ggridges_0.5.3
## [69] Rcpp_1.0.6 sparseMatrixStats_1.2.0
## [71] zlibbioc_1.36.0 RCurl_1.98-1.2
## [73] ps_1.5.0 rpart_4.1-15
## [75] deldir_0.2-9 viridis_0.5.1
## [77] pbapply_1.4-3 zoo_1.8-8
## [79] haven_2.3.1 ggrepel_0.9.1
## [81] cluster_2.1.1 fs_1.5.0
## [83] magrittr_2.0.1 RSpectra_0.16-0
## [85] scattermore_0.7 openxlsx_4.2.3
## [87] lmtest_0.9-38 reprex_1.0.0
## [89] RANN_2.6.1 fitdistrplus_1.1-3
## [91] hms_1.0.0 mime_0.10
## [93] evaluate_0.14 xtable_1.8-4
## [95] readxl_1.3.1 compiler_4.0.3
## [97] KernSmooth_2.23-18 crayon_1.4.1
## [99] htmltools_0.5.1.1 mgcv_1.8-33
## [101] later_1.1.0.1 lubridate_1.7.10
## [103] DBI_1.1.1 formatR_1.7
## [105] tweenr_1.0.1 ExperimentHub_1.16.0
## [107] dbplyr_2.1.0 MASS_7.3-53.1
## [109] rappdirs_0.3.3 Matrix_1.3-2
## [111] cli_2.3.0 rbibutils_2.0
## [113] metap_1.1 igraph_1.2.6
## [115] pkgconfig_2.0.3 plotly_4.9.3
## [117] scuttle_1.0.0 xml2_1.3.2
## [119] vipor_0.4.5 XVector_0.30.0
## [121] rvest_1.0.0 digest_0.6.27
## [123] sctransform_0.3.2 RcppAnnoy_0.0.18
## [125] spatstat.data_2.0-0 rmarkdown_2.6
## [127] cellranger_1.1.0 leiden_0.3.7
## [129] uwot_0.1.10 DelayedMatrixStats_1.12.0
## [131] curl_4.3 shiny_1.6.0
## [133] lifecycle_1.0.0 nlme_3.1-152
## [135] jsonlite_1.7.2 BiocNeighbors_1.8.0
## [137] fansi_0.4.2 limma_3.46.0
## [139] viridisLite_0.3.0 pillar_1.4.7
## [141] lattice_0.20-41 fastmap_1.1.0
## [143] httr_1.4.2 survival_3.2-7
## [145] interactiveDisplayBase_1.28.0 glue_1.4.2
## [147] zip_2.1.1 spatstat_1.64-1
## [149] png_0.1-7 BiocVersion_3.12.0
## [151] bit_4.0.4 ggforce_0.3.2
## [153] stringi_1.5.3 blob_1.2.1
## [155] BiocSingular_1.6.0 AnnotationHub_2.22.0
## [157] memoise_2.0.0 irlba_2.3.3
## [159] future.apply_1.7.0writeLines(capture.output(sessionInfo()), "./scripts/9_Subcluster_Macrophages/9_Subcluster_Macrophages.sessionInfo.txt")